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Science  and  Engineering 


This  thesis  develops  a  program  which  will  merge  or 
overlay  imagery  and  terrain  elevation  data  and  create  a 
synthetic  3-D  perspective  view  of  the  ocean  bottom.  The 
observer  may  position  himself  at  various  locations  and  see 
the  terrain  from  different  viewpoints.  The  elevation  data 
is  grouped  into  triangular  panels  and  the  color 
information  is  averaged  from  the  imagery  data  file.  The 
entire  panel  is  assigned  a  single  color  equal  to  the 
average.  These  panels  are  then  projected  onto  an  image 
plane  by  using  a  3D  to  2D  perspective  transformation. 
Hidden  surfaces  are  removed  by  a  , "painters”  algorithm 
which  relies  on  sorting  the  panels  based  on  distance  from 
the  observer. 
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I.  INTRODUCTION  fl 

The  goal  of  this  thesis  is  to  develop  the  software 
required  to  display  a  3-D  perspective  view  of  combined 
sonar  imagery  and  bathymetric  data.  The  result  is  a 
synthesized  image  representing  a  3-D  perspective  view  of 
the  terrain  from  a  given  observer  location.  I 

The  issues  involved  include:  first,  how  can  the  sonar 
imagery  data  and  the  bathymetric  data  be  combined;  second, 
what  is  the  methodology  of  creating  the  3-D  perspective 
view  from  all  observer  locations;  third,  how  fast  can  each 
view  be  generated  using  real  data  and  fourth,  how  is  the 
picture  quality  affected  by  the  resolution  of  the  data. 

Conventional  methods  of  displaying  elevation  data  are 
with  contour  lines  or  3-D  grid  line  drawings.  Some 
applications  merge  the  contour  lines  with  color  data  and 
this  aids  in  the  overall  comprehension  of  the  data. 

Recently  the  speed  of  the  digital  computer  has  been  called 
upon  to  generate  3-D  views  based  on  elevation  with 
shading.  The  shading  value  may  be  based  on  the  elevation, 
or  may  be  based  on  other  information  available.  These 
displays  are  an  improvement  over  the  simple  contour  plots, 
as  well  as  the  3-D  grid  line  drawings.  The  advantage  of 
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any  3-D  perspective  display  method  is  the  ability  to 
"see"  the  terrain  in  a  fashion  which  appears  normal  to 
the  observer. 

Applications  of  this  method  range  from  low  cost  flight 
simulators  for  manned  or  unmanned  vehicles  to  real  time 
displays  of  multiple  source  information  in  a  fashion  which 
results  in  a  greater  level  of  comprehension  or 
interpretation.  For  example,  a  combat  aircraft  pilot  must 
be  able  to  assimilate  vast  amounts  of  data.  The  pilot  must 
search  multiple  electronic  displays  in  order  to  put 
together  a  picture  of  the  environment  in  his  mind.  The 
pilot  must  answer  several  questions  simultaneously:  Where 
am  I?,  Where  are  my  friends?.  Where  is  the  enemy?,  What 
is  the  condition  of  my  aircraft?  Simplifying  the 
presentation  of  this  information  in  a  simulated 
perspective  view  of  the  surrounding  environment  will 
improve  the  response  of  the  pilot.  [Ref .  1,  p.  64]  The 
ability  to  combine  the  various  data  types  into  a  combined 
display  is  a  form  of  "sensor  fusion".  In  the  aircraft 
example  the  data  may  be  video  from  external  video  cameras 
and  radar  or  infrared  sensor  returns. 

Generally,  sensor  fusion  refers  to  the  ability  to 
merge  layers  of  information  from  various  sources  into  a 
new  form.  The  data  may  be  imagery  data  from  Landsat  or 
aircraft,  or  any  other  layer  of  data  within  the  same 
geographical  area  such  as  magnetic  anomaly  or  infrared 
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response  data.  These  various  types  of  data  must  be  in  a 
form  which  can  be  represented  by  a  specific  crdcr  or  gray 
shade.  Adding  color  or  shading  from  the  other  data  source 
allows  the  correlation  between  it  and  the  terrain 
elevation  data  to  become  apparent.  A  specific  area  of 
interest  includes  shipboard  systems  where  numerous  sensors 
are  gathering  data  and  the  operator  selectively  views  the 
output  from  a  single  sensor  or  a  combination  of  sensors. 
If  the  outputs  from  multiple  sensors  are  combined  into  a 
3-D  perspective  view  of  the  surrounding  environment,  the 
response  time  of  the  operator  will  improve.  By  utilizing 
the  3-D  techniques,  a  view  forward  of  the  ship  could  be 
displayed  on  a  CRT  screen.  Combining  radar  elevation  data 
with  an  image  file  database  would  permit  the  formation  of 
a  3-D  perspective  view  of  the  radar  image.  This  same 
technique  used  with  the  sonar  data  obtained  by  the  ship's 
sonar  systems  and  imagery  data  files  of  the  harbor  channel 
would  allow  a  realistic  view  of  the  channel  and  permit 
more  efficient  usage  of  the  available  sonar  systems. 

The  software  developed  in  this  thesis  overlay  side¬ 
looking  sonar  imagery  data  of  the  ocean  bottom  onto  the 
bathymetric  data  of  the  same  geographic  area.  This  is  an 
extension  of  work  by  L.  Coleman  [Ref.  2].  Coleman's  work 
concentrated  in  generating  a  3-D  perspective  view  of  land 
terrain.  Several  items  added  in  this  work  include  the 
ability  to  approach  the  area  from  any  heading  (observer 
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locations  within  the  area  boundaries  are  permitted) , 
reduction  of  the  time  to  create  one  image,  and  drawing  of 
the  image  is  included  in  the  main  routine.  In  addition  to 
the  work  by  Coleman,  a  separate  project  by  McGhee,  Zyda, 
Smith  and  Streyle  incorporates  many  of  the  same  techniques 
[Ref.  3].  The  imagery  data  is  in  the  form  of  a  512  x  512 
pixel  image  with  256  possible  gray  levels.  The  bathymetric 
data  is  gridded  and  consists  of  depth  values  on  a  latitude 
longitude  grid.  The  bathymetric  data,  called  the  elevation 
data,  is  segmented  into  triangular  panels.  The  boundaries 
of  each  panel  overlay  the  image  data,  image  pixels  inside 
the  panel  are  averaged  and  this  average  gray  level  is 
assigned  to  the  entire  panel .  To  generate  the  perspective 
view,  each  elevation  point  is  projected  onto  an  image 
plane  located  at  the  observer  location,  and  oriented  with 
respect  to  the  observer  course.  This  projection  is 
accomplished  through  a  3D-2D  perspective  projection 
transformation  which  maps  the  elevation  point  coordinates 
to  image  plane  coordinates.  The  synthesized  image  is 
generated  on  a  monitor  and  approximates  the  view  from  the 
observer  position  with  a  38.6  degree  field  of  view.  The 
synthesized  view  is  directly  affected  by  the  resolution  of 
the  input  data.  The  elevation  file  resolution  affects  both 
the  elevation  drawing  and  the  shading  of  the  image. 
Shading  is  affected  due  to  the  averaging  of  pixel  shades 
within  each  triangle.  The  drawing  is  affected  by  the  low 


sampling  rate  of  the  terrain;  only  the  larger  features 
will  be  present  in  the  displayed  image. 


The  software  is  implemented  in  FORTRAN  on  a  DEC  Micro- 
Vax  GPX  II.  This  system  is  capable  of  1024  x  864  pixel 
resolution,  with  256  colors  or  256  gray  levels.  A  graphics 
software  package  from  DEC,  MicroVMS  Workstation  Graphics 
Software  version  3.0,  is  used  to  draw  and  fill  each  panel 
on  the  monitor.  The  software  may  be  run  on  other  systems, 
with  changes  required  in  the  screen  plotting  routines  and 
the  number  of  shades  or  color  selection.  Also  the  size  of 
the  input  data  files  may  be  limited  on  other  systems  which 
do  not  have  sufficient  memory  to  store  the  input  data  in 
arrays.  Results  of  this  work  show  that  the  performance  is 
limited  by  the  plotting  speed  of  the  Micro-Vax.  The 
calculation  of  the  perspective  image  is  approximately  11% 
of  the  total  time  required  to  complete  the  process.  The 
remainder  of  time  is  used  to  draw  the  image  on  the  screen. 

Using  this  software,  an  observer  may  select  a  position 
and  heading  and  "see”  the  ocean  bottom  terrain  in  a  3-D 
perspective  form.  The  observer  heading  is  not  limited,  the 
position  may  be  within  the  boundaries  of  the  area 
selected,  and  the  view  is  generated  directly  on  the 
screen.  Several  options  are  available  including  a 
magnification  factor,  elevation  scale  exaggeration,  saving 
the  image  to  disk  and  the  ability  to  respecify  a  new 
observer  location.  The  speed  of  image  generation  is 
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directly  affected  by  the  resolution  or  size  of  the  input 
data  files,  using  61  x  61  elevation  points,  an  image  is 
generated  in  40  seconds. 

Chapter  Two  will  cover  the  basic  image  geometry  in 
order  to  develop  the  perspective  projection  transformation 
equation,  and  the  orientation  of  the  image  plane.  Chapter 
Three  details  the  data  file  formats,  resolution  and  how 
they  are  prepared  for  use  in  the  main  routine.  Chapter 
Four  covers  the  implementation  of  the  algorithms  in  the 
program  and  any  alternatives  which  were  investigated. 
Chapter  Five  details  program  operation  and  Chapter  Six 
presents  conclusions  and  several  recommendations  for 
improvement  or  further  study. 


II.  IMAGE  FORMATION  USING  TRANSFORMATION  GEOMETRY 


Inherent  in  this  project  is  the  ability  to  take  the 
two  dimensional  gridded  elevation  data  which  represents  a 
three  dimensional  object  and  display  a  perspective  view  on 
a  two  dimensional  monitor.  Two  basic  approaches,  namely, 
parallel  or  perspective  projection  can  accomplish  this. 

A  parallel  projection  is  one  where  the  points  of  the 
object  are  projected  onto  the  image  plane  by  parallel 
lines,  that  is  the  projection  lines  do  not  converge.  This 
method  has  the  advantage  of  maintaining  relative 
dimensions  of  the  object  and  is  useful  for  obtaining 
measurements  [Ref.  3,  p.  52]. 

In  the  perspective  projection  all  points  are  projected 
onto  the  image  plane  through  a  reference  point  which  will 
be  called  the  focal  point  [Ref.  4,  p.  133].  In  this 
method,  the  relative  dimensions  of  the  object  are  not 
preserved,  but  the  image  appears  more  realistic.  Figure 
2-1  illustrates  the  two  methods.  Because  the  goal  is 
feature  recognition  and  interpretation  through  screen 
display  and  not  precise  measurement,  the  perspective 
projection  is  used. 
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Figure  2-1  Parallel  and  Perspective  Projections 

[Ref.  3,  p.  54) 


The  perspective  projection  may  be  generalized  as  a 
reference  3-D  to  2-D  transformation  and  is  illustrated  in 
matrix  notation: 

A  =  [  M  ]  B  (2.1) 

where:  A-represents  a  vector  in  the  image  3D 
coordinate  system 

B-represents  a  vector  in  the  object  3D 
coordinate  system 

M-represents  the  transformation  matrix 

A.  COORDINATE  SYSTEMS 

Equation  2.1  refers  to  the  image  3-D  coordinate  system 
and  the  object  3-D  coordinate  system.  The  object  is 
located  in  a  right-handed  3-D  cartesian  coordinate  system 
where  each  object  point  is  described  in  terms  of  (X,Y,Z) 
geo-rectangular  coordinates.  The  center  of  the  earth  is 
located  at  (0,0,0),  the  X-axis  points  to  the  intersection 
of  0  degrees  latitude  and  0  degrees  longitude,  the  Y- 
axis  points  to  the  intersection  of  0  degrees  latitude  and 
90  degrees  east  longitude  and  the  Z-axis  points  to  true 
north.  Figure  2-2  illustrates  this  coordinate  system. 
There  is  a  direct  relationship  between  the  latitude, 
longitude,  height  with  respect  to  sea  level  and  the 
(X , Y , Z )  object  coordinates.  The  input  data  is  supplied  on 
a  latitude,  longitude  grid  and  is  converted  to  (X,Y,Z) 
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coordinates  for  use  in  the  transformation  equations.  The 
image  coordinate  system  is  also  a  right-handed  3-D 
cartesian  system.  The  focal  point  is  located  at 
(xO,yO,F).  Figure  2-3  illustrates  the  image  coordinate 
system. 

B.  3-D  PERSPECTIVE  TRANSFORMATION 

The  elements  of  the  [M]  matrix  in  equation  2.1  are  the 
cosines  of  the  spatial  angles  that  each  of  the  image 
system  axis  x,y,z  make  with  each  of  the  object  system  axis 
X, Y , Z .  This  substitution  results  in  the  following  specific 
form  for  the  matrix  [M] .  [Ref.  4,  p.  139] 
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cosxX  cosxY  cosxZ 

M11  M12  m13 

M  = 

cosyX  cosyY  cosyZ 

h21  m22  m23 

coszX  coszY  coszZ 

M3 1  M32  M33 

»  - 

By  convention  the  [M]  matrix  transforms  from  the 
object  system  to  the  image  system  [Ref.  4,  p.  139].  Using 
this  matrix,  every  object  point  may  be  converted  from  the 
3-D  XYZ  coordinates  to  the  xyz  image  coordinates. 

At  this  point  the  relation  for  a  single  object  point 
is  developed.  Once  defined,  the  relation  may  then  be 
applied  to  all  the  vectors  from  the  object  to  the  focal 
point.  Figure  2-4  illustrates  the  object  and  image  space 
coordinate  systems.  Define  the  vector  <FA>  as  a  vector 
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Figure  2-4  Object  and  Image  Space  Systems 


from  the  focal  point  (F)  to  a  point  on  the  image  plane  (A) 
and  the  vector  <FB>  as  a  vector  from  the  focal  point  (F) 
to  an  object  point  (B)  .  Note  that  the  image  vector  <FA> 
is  collinear  with  the  object  vector  <FB>.  [Ref.  4,  p. 
141] 

FA  =  k  FB 

where:  FA  -  image  vector 

FB  -  object  vector  (2.3) 

k  -  constant  of  proportionality 
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Using  the  image  space  coordinates  (x,y,z)  to  describe  <FA> 
and  object  space  coordinates  (X,Y,Z)  to  describe  <FB> 


I 


gives: 


xA  -  x0 

'  ;B 

FA  - 

ya  -£Yb 

FB  « 

zB 

ZB 

'  zF 
-  zF  J 

(2.4) 


Then  using  equation  2.1  to  express  <FA>  in  the  xyz  image 
system  yields: 

FA  -  k  [  M  J  FB  (2.5) 

or  after  substitution: 


XA  '  x0 

'  xB  -  xF  1 

Ya  -  Y0 

'  k  [  M  ] 

yb  “  yf 

-£ 

1 

ZB  -  Zp 

Equation  2.6  is  a  form  of  the  collinearity  equation,  it 
forms  the  basic  relationship  that  the  focal  point,  the 
image  point  and  the  object  point  are  in  a  straight  line. 
The  values  of  xO  and  yO  allow  for  the  observer  to  be 
slightly  misaligned  with  the  image  plane  z-axis  and 
generally  are  set  equal  to  zero.  This  is  the  fundamental 
relationship  used  in  this  project  to  create  the 
perspective  views  on  the  monitor  screen.  [Ref.  4,  p.  141] 
This  transformation  will  be  called  upon  after  the 
elevation  points  have  been  grouped  into  triangles  and  the 
observer  location  has  been  entered.  Each  elevation  point 
in  front  of  the  observer  will  be  transformed  to  image 
plane  coordinates  for  subsequent  drawing  on  the  2-D 
screen. 
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III.  DATA  FILE  FORMAT  FOR  DIFFERENT  SENSOR  IMAGES 

Input  to  the  routines  consist  of  two  data  files,  the 
bathymetric  data  and  the  image  data.  Each  of  these  is 
described  below. 

A.  BATHYMETRIC  DATA 

The  bathymetric  data  was  obtained  from  the  Defense 
Mapping  Agency/National  Geophysical  Data  Center.  The 
complete  data  base  should  be  referenced  as  "Digital 
Bathymetric  Data  Base-Unclassified"  or  (DBDBU) .  The  area 
which  was  used  here  extends  from  50  degrees  north 
latitude,  140  west  longitude  to  32  degrees  north,  120  west 
on  a  5  minute  by  5  minute  grid.  Each  value  describes  the 
depth  in  meters  below  sea  level  at  a  given  coordinate 
location.  Values  which  lie  above  sea  level  are  assigned 
-10  to  avoid  ambiguity. 

Extracting  the  area  of  interest  is  performed  by  a 
utility  program,  CROPELEV.  FOR.  The  user  enters  the 
latitude  and  longitude  of  the  southwest  corner  of  the  area 
of  interest  and  the  number  of  rows  and  columns  desired. 
Each  row  and  column  is  5  minutes  of  latitude  and  longitude 
respectively. 


After  the  area  of  interest  is  selected,  the  data  is 
extracted  from  the  original  file  and  placed  in  row- 
column  format  where  the  rows  correspond  to  latitude  and 
the  columns  represent  the  longitude.  Row  numbers  increase 
from  south  to  north  and  column  numbers  increase  from  west 
to  east.  The  first  record  of  the  file  is  a  header 
containing  the  latitude  and  longitude  (deg,  min,  sec)  of 
the  southwest  corner,  the  latitude  and  longitude 
resolution  in  seconds  and  the  number  of  rows  and  columns 
of  data.  Figure  3-1  illustrates  this  file  format.  This 
file  is  stored  on  disk  and  loaded  into  the  variable 
RELEV (*,*)  at  runtime. 


B.  IMAGERY  DATA 

The  sonar  imagery  data  was  obtained  from  the  U.S. 
Geological  Survey,  Department  of  the  Interior  "Atlas  of 
the  Exclusive  Economic  Zone,  Western  Conterminous  United 
States",  [Ref.  5].  This  atlas  consists  of  36  two  degree 
mosaic  images  at  a  scale  of  1:500,000.  Coverage  extends 
from  49  degrees  north,  13  0  west  to  3  0  degrees  north,  117 
west.  These  images  were  obtained  using  a  unique  side-scan 
sonar  system  called  GLORIA  (Geological  Long-Range  Inclined 
Asdic).  [Ref.  5,  p.  2]  This  sonar  system  allows  mapping 
of  large  areas  of  ocean  bottom  on  a  single  pass  of  the 
ship.  Figure  3-2  lists  several  characteristics  of  the 
GLORIA  system. 
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Figure  3-1  Bathymetric  Data  File  Format 
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Gloria  Side  Scan  Sonar  Specifications 


Size  .  7.75  x  .66  meters 

Weight  .  2  tons 

Scan  .  30  km/side  =  60  km 

swathwidth  at  10  kts. 

Power  .  10  kw/side 

Beamwidth  .  2.5  deg  azimuth 

10  deg  vertical 

Resolution  .  30  x  218  meters/pixel 

at  5000  meters  depth 

Frequency  .  6.5  kHz 


Figure  3-2  Some  Characteristics  of  the 
Gloria  II  Side-Scan  Sonar 
[Ref.  6,  p.  7] 


Each  mosaic  image  is  a  half-tone  black  and  white 
print  of  the  acoustic  reflectance  of  the  sea  floor  with 
white  representing  the  highest  reflectance  and  black  the 
lowest  reflectance.  As  seen  in  Figure  3-2,  the  resolution 
of  each  pixel  is  distorted  (30  meters  by  218  meters)  due 
to  the  ships  motion  along  a  track  perpendicular  to  the 
scanning  direction.  The  smaller  value  (30  meters)  is  the 
resolution  in  the  cross-track  direction  while  the  larger 
value  (218  meters)  is  in  the  along  track  direction.  Prior 
to  forming  the  mosaic  images,  aspect  ratio  distortion  is 
removed  by  correcting  for  geometric  and  radiometric 


17 


distortions  in  the  raw  data.  The  final  resolution,  after 
correction,  used  to  create  the  mosaics  is  about  50  meters 
x  50  meters.  This  relatively  low  resolution  image  can 
show  only  large  scale  features  but  may  point  out  areas 
which  deserve  additional  attention.  Also  note  that  this 
resolution  is  several  orders  of  magnitude  better  than  the 
available  bathymetric  data.  [Ref.  6,  p.  3] 

Digital  sonar  imagery  data  was  not  available  for  use 
in  this  thesis.  In  order  to  determine  the  effectiveness  of 
this  imaging  method,  the  digital  image  data  was  created  by 
locally  digitizing  the  mosaics  contained  in  the  atlas. 
This  provided  digital  imagery  data  in  the  following  form, 
512  x  512  pixels  with  each  pixel  using  one  byte  of  storage 
resulting  in  256  shades  of  gray.  This  data  is  stored  on 
disk  in  direct  access  form  as  512  fixed  length  records 
with  512  bytes  per  record.  Record  one  is  the  northern  end 
of  the  image  and  record  512  is  the  most  southern  end  of 
the  image.  Figure  3-3  illustrates  the  orientation  of  the 
image  file.  The  disk  file  is  loaded  into  the  variable 
IM\GE ( * , *)  at  runtime. 

Ideally,  the  digital  imagery  data  would  be  directly 
available,  either  from  the  GLORIA  side-scan  sonar  or  from 
other  sonar  systems.  If  this  were  the  case,  the  data  would 
require  processing  to  remove  the  different  sources  of 
error  and  to  convert  it  from  its  native  form  into  the 
form  required  in  this  project.  Processing  techniques  for 
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GLORIA  digital  data  is  discussed  in  Ref.  7,  "Processing 
Techniques  for  Digital  Sonar  Images  from  GLORIA"  by  Pat  S. 
Chavez,  Jr. 

The  current  limitation  for  input  data  is  70  rows  by  70 
columns  for  the  elevation  data  and  1024  rows  by  1024 
columns  for  the  image  data.  These  dimensions  may  be 
increased  to  accommodate  higher  resolution  data;  the 
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tradeoff  is  speed  of  execution.  The  number  of  triangular 
panels  is  a  function  of  the  elevation  data  size, 
specifically: 


#panels  =  ( ( #rows  -  1)*  2)  *  (#cols  -  1)  (3.1) 


An  elevation  file  with  resolution  20  x  20  has  722 
panels;  increasing  the  resolution  to  100  x  100  results  in 
19,602  panels.  Each  panel  is  projected  point  by  point  then 
sorted  and  drawn  for  every  image  view  constructed.  The 
increase  in  computation  time  follows  directly.  In  addition 
to  the  elevation  file  resolution,  the  image  file 
resolution  may  be  increased.  The  cost  in  computation  is 
not  severe  in  this  case  because  the  image  file  is  accessed 
only  once  when  the  image  pixels  are  averaged  to  calculate 
the  color  or  gray  level  for  the  individual  panels.  This  is 
done  prior  to  displaying  the  first  image  and  is  not 
repeated  unless  the  program  is  exited  and  restarted. 
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IV.  ALGORITHMS  FOR  COMBINED  3-D  PERSPECTIVE  DISPLAY 


The  main  routine  begins  by  reading  two  input  data 
files  into  the  arrays  RELEV(*,*)  and  IMAGE(*,*).  Figure 
4-1  illustrates  the  steps  involved  in  displaying  the  image 
on  the  monitor. 

The  procedure  begins  by  forming  the  triangles  in  the 
grid  of  the  elevation  data,  averaging  the  image  pixels 
inside  each  triangle  and  assigning  the  average  gray 
intensity  to  the  panel  (one  triangle) .  The  observer 
location  is  read,  then  the  image  plane  coordinate  system 
is  constructed.  Based  on  the  orientation  of  the  image 
plane  the  [M]  matrix  is  calculated  and  used  to  map  the 
elevation  points  to  the  image  plane.  A  2-D  to  2-D  affine 
transform  is  used  to  convert  the  points  from  image  plane 
coordinates  to  the  screen  coordinates.  To  display  the 
perspective  view  on  a  flat  screen,  triangles  which  are 
hidden  behind  foreground  objects  must  be  removed.  This 
hidden  surface  removal  is  performed  by  calculating  the 
distance  from  the  panel  to  the  observer  and  drawing  the 
panels  in  sequence  from  the  farthest  to  the  nearest.  This 
chapter  will  cover  these  areas  in  detail  and  also  discuss 
alternatives  which  were  considered,  but  not  implemented. 
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Begin 


Load  input  data  file 

Read  in  elevation  (bathymetric)  data 
Read  in  imagery  data 

Construct  triangles  in  elevation  data  grid 
Average  the  gray  shade  in  each  triangle 

Get  observer  location 

Form  the  image  plane  vectors 

Calculate  the  elements  of  the  [M]  matrix 

Transform  to  the  image  plane 

Determine  which  panels  are  visible 

Convert  to  screen  coordinates 

Remove  hidden  surfaces 

Plot  results  on  screen 

End 

Figure  4-1  Basic  Program  Flow 


A.  POLYGON  FORMATION  AND  SHADING 

Solid  objects  may  be  represented  in  numerous  ways. 
Some  objects  lend  themselves  to  being  described  in  terms 
of  a  number  of  planes  or  surfaces.  A  cube,  for  example  may 
be  precisely  defined  with  six  planes  [Ref.  8,  p.  189].  As 
the  object  or  scene  becomes  more  complex,  the  number  of 
surfaces  required  to  accurately  describe  it  increases. 
This  method  of  representing  a  3-D  surface  by  plane 
surfaces  lends  itself  particularly  well  to  this 
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application.  It  allows  for  easy  calculation  of  the  color 
or  gray  level  by  averaging  of  pixels  contained  within  the 
panel,  easy  transformation  to  the  image  plane,  and  easy 
use  of  graphics  hardware  to  execute  the  polygon  fill 
operation.  This  last  feature  is  most  important;  earlier 
work  in  this  area  pointed  out  the  amount  of  time  required 
to  draw  based  on  a  point  by  point  method  [Ref.  2,  p.  56]. 
Using  the  graphics  hardware  polygon  fill  capability 
alleviates  this  problem. 

The  input  elevation  data  is  supplied  in  gridded  form 
and  is  easily  partitioned  into  squares  and  ultimately 
into  triangles.  Figure  4-2  illustrates  the  partitioning  of 
the  elevation  data  into  triangles.  The  procedure  for 
selecting  the  gridsquares  and  forming  the  triangles  is 
given  as  psuedocode  in  Figure  4-3.  The  column  coordinates 
of  each  point  are  stored  in  IA(*)  and  the  row  coordinates 
are  stored  in  the  JA(*)  array.  Starting  in  the  southwest 
corner,  the  program  selects  four  elevation  points  which 
form  a  gridsquare.  These  four  points  are  called  node_a, 
node_b,  node_c  and  node_d.  The  gridsquare  is  divided  into 
two  triangles  by  calculating  the  equation  of  the  line 
which  divides  it.  This  dividing  line  is  completely 
described  by  the  slope  and  y-intercept.  The  column 
boundaries  are  defined  by  IA(node_a)  and  IA(node_b) ;  the 
row  boundaries  are  defined  by  JA(node_a)  and  JA(node_c) . 
For  each  incremental  step  from  the  eastern  (right  hand) 
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North 

I 


View  from  above  looking  down  on  the  terrain. 

-Terrain  elevation  points  are  connected 
to  form  triangular  polygons  with  common 
edges . 


Figure  4-2  Polygonal  Terrain  Construction 

[Ref.  3,  p.  35] 
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BEGIN 


DO  select  a  gridsquare 

(♦select  four  nodes  which  make  one  gridsquare*) 

READ  node_a 

READ  node_b 

READ  node_c 

READ  node_d 

(*calc  slope  and  y-intercept  of  line  forming 
the  triangle*) 

slope  =  rise  /  run 

y-intercept  =  y  -  slope  *  (x) 

IX)  M  =  IA(node_b)  to  IA(node_a)  step  -1 

IY=( slope)  *  M  +  y-intercept 

DO  L  *  JA(node_b)  to  IY  step  -1 

(♦average  image  pixels  in  lower  triangle*) 

END  DO 

DO  L  =  IY  to  JA(node_c)  step  -1 

(♦average  image  pixels  in  upper  triangle*) 
END  DO 
END  DO 

PL_AR(*,*)  =  nodes  of  triangle,  average  gray  shade 
END  DO 
END 

Figure  4-3  Polygon  Formation  Psuedocode 
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boundary,  IA(node_b)  to  the  western  (left  hand)  boundary, 
IA(node_a)  ,  a  value  of  y  equal  to  IY  is  calculated  based 
on  the  equation  of  the  dividing  line.  The  image  pixels 
above  and  below  this  line  are  averaged  separately  for  the 
upper  and  lower  triangles.  The  process  is  repeated  for 
each  incremental  step  in  the  columns  from  right  to  left. 
This  scanning  and  averaging  pattern  is  illustrated  in 
Figure  4-4.  Note  that  the  scan  pattern  is  from  bottom  to 
top,  right  to  left  in  each  triangle.  [Ref.  2,  p.  39]  The 
average  shade  of  the  triangle  along  with  the  panel  number 
and  vertices  are  stored  in  array  PL_AR(*,*).  The  data 
structure  is  shown  in  Figure  4-5.  Notice  that  PL_AR(*,5) 
generally  is  not  used,  it  will  be  used  to  store  the 
maximum  distance  of  the  plane  from  the  observer. 

B.  IMAGE  PLANE  FORMATION,  PERSPECTIVE  VIEW  CALCULATION 

After  the  shading  for  each  panel  is  calculated,  the 
observer's  location  is  entered.  This  consists  of  latitude 
and  longitude  (deg,  min,  sec),  height  with  respect  to  sea 
level  and  heading.  The  latitude,  longitude  and  height  data 
is  converted  to  object  space  coordinates  X1,Y1,Z1  while 
the  heading  is  used  to  select  a  second  point  in  front  of 
the  observer.  This  second  point  X2,Y2,Z2  forms  the  line  of 
sight  (LOS)  vector  for  the  observer.  The  negative  LOS 
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Figure  4-5  Layout  of  PL_AR  Array 


vector  is  the  z-axis  of  the  image  coordinate  system  as 
seen  in  Figure  2-3.  This  is  formulated  as  shown: 

ZVECX  -  XI  -  X2 

ZVECy  =  Y1  -  Y2  (4.1) 

ZVECZ  =  Z1  -  Z2 

ZVEC  =  ZVECX  x  +  ZVECy  y  ♦  ZVECZ  z 
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The  image  plane  y-axis  is  constructed  by  creating  a  vector 
which  points  upward  from  the  center  of  the  earth  to  the 
observer  location.  The  (X,Y,Z)  coordinates  of  the  observer 
location  (X1,Y1,Z1)  form  this  vector.  The  y-axis  is 
described  by: 


(4.2) 


YVECX  =  XI 
YVECy  »=  Y1 
YVECZ  =  21 

YVEC  =  YVECX  x  +  YVECy  y  +  YVECZ  z 


The  image  coordinate  system  is  a  right-hand  system. 
Therefore,  the  cross  product  of  the  y-axis  and  z-axis 
forms  the  x-axis.  The  x-axis  of  the  image  plane  is 
constructed  as  follows: 


(4.3) 


XVECX  «  ((YVECy  x  ZVECZ)  -  (YVECZ  x  ZVECy)) 
XVECy  =  ((YVECZ  x  ZVECX)  -  ( YVECX  x  ZVECZ)) 
XVECZ  =  ((YVECX  x  ZVECy)  -  (YVECy  x  ZVECX)) 

XVEC  =  XVECX  x  +  XVECy  y  +  XVECZ  z 


The  [M]  matrix  is  calculated  using  the  relationship 
developed  in  Chapter  Two,  equation  2.2.  Recall  that  each 
element  of  the  [M]  matrix  is  the  cosine  of  the  spatial 
angle  between  the  axis  of  the  image  coordinate  system  and 
the  object  coordinate  system  axis.  The  object  space 
coordinate  axis  may  be  represented  as  unit  vectors: 


XAXIS  **  lx  +  0y  +  0z 
YAXIS  «  0x  +  ly  +  0z 
ZAXIS  =  0x  +  0y  +  lz 


(4.4) 
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The  cosine  of  the  spatial  angle  between  two  vectors  is 
determined  from  the  dot  product  of  the  two  vectors. 

i 

A  •  B  *  |  A  |  |  B  |  cos  ©AB  (4.5) 


The  dot  product  may  be  expanded  as  shown: 


cos  ©ab 


A  ‘  B  . 
IAI  IB  I 


cos  ©AB 


AgBg  *  Ay By  +  A2B2 

IAI  IB  I 


(4.6) 


Substituting  the  image  space  x-axis  vector  (xvec)  and  the 
object  space  X-axis  vector  (XAXIS)  for  [A]  and  [B]  allows 


Mil  to  be  calculated: 


Mu*  cos  exX 


XVECy -XAXISy  +  XVECy "XAXISy  +  XVEC2"XAXIS2 
IXVECI  IXAXIS  I 


(XVECX -1)  +  (XVECy-0)  +  (XVEC2-0) 
IXVECI  •  1 


(4.7) 


XVECy 

Hit  ■  ■■  ■  — 

11  IXVECI 

The  other  elements  of  the  [M]  matrix  may  be  found  in  a 
similar  fashion.  This  results  in: 


M11  “ 

XVECy 
IXVEC  | 

m12  * 

XVECy 
IXVEC  1 

m13  * 

XVEC2 
IXVEC  1 

M21  " 

YVECy 

IYVECI 

M22  * 

YVECy 
IYVEC  1 

M23  “ 

yvec2 

IYVEC  1 

(4.8) 

m3i  - 

ZVECy 
IZVEC  1 

m32  * 

ZVECy 
IZVEC  1 

m33  * 

ZVECZ 
IZVEC  1 
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This  method  of  calculating  the  image  coordinate  axis 
and  the  corresponding  [M]  matrix  elements  results  in  a 
viewing  plane  which  is  oriented  vertically  and  faces  in 
the  direction  of  the  observer  course. 

An  alternative  to  this  method  is  that,  instead  of 
specifying  the  course  to  see  the  perspective  view,  one  can 
have  the  image  plane  always  face  the  object.  This  would 
simulate  rotating  the  object  while  keeping  the  observer 
location  fixed.  The  decision  to  use  the  first  method  is 
based  on  the  feeling  that  it  results  in  a  more  natural 
view  for  the  observer. 

C.  TRANSFORM  TO  THE  IMAGE  PLANE 

After  the  [M]  matrix  is  constructed  the  elevation 
points  are  projected  onto  the  image  plane.  Points  which 
lie  behind  the  image  plane  or  behind  the  observer  are  not 
projected.  The  method  used  to  select  which  points  to 
project  and  the  projection  equations  are  presented  next. 

Determining  which  elevation  points  are  located  behind 
the  image  plane  is  done  by  finding  the  cosine  of  the 
angle  between  a  vector  formed  by  the  negative  z-axis  of 
the  image  plane  and  a  vector  drawn  from  X2,Y2,Z2  to  the 
object  point.  Solving  for  the  cosine  of  the  angle  and  not 
the  angle  itself  eliminates  the  need  to  use  trigonometric 
functions,  and  provides  faster  execution.  Psuedocode  for 
this  procedure  is  given  in  Figure  4-6. 
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BEGIN 

DO  IR  =  1  to  Number  of  elevation  pts. 

(*get  XYZ  coordinates  of  point*) 

X  =  XYZ (IR, 1) 

Y  =  XYZ (IR, 2) 

Z  =  XYZ (IR, 3) 

(*form  the  observer  LOS  vector*) 

OBS_VECX  =  X2  -  XI 
OBS_VECY  =  Y2  -  Y1 
OBS_VECZ  =  Z2  -  Z1 

(♦form  the  object  vector*) 

OBJ_VECX  =  X  -  X2 
OBJVECY  =  Y  -  Y2 
OBJ_VECZ  =  Z  -  Z2 

(*find  the  cosine  of  the  angle  between  the 
OBS_VEC  and  OBJ_VEC* ) 

COS_THETA  «  (dot  product  of  OBS_VEC,  OBJ_VEC)/ 
(Mag (OBS_VEC)  x  Mag (OBJ_VEC) ) 

IF  ( COS_THETA  >  0)  THEN 

Project  onto  image  plane 
Calculate  depth  of  panel 

ELSE 

Mark  as  a  non-visible  point 
END  IF 


END  DO 

END 


* 


Figure  4-6  Psuedocode  Determine  Which  Panels  to  Project 
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The  procedure  begins  by  forming  two  vectors,  the 
observer  LOS  vector  and  the  observer  to  object  vector 
(object  vector) .  Ideally,  the  origin  of  the  image  plane 
coordinate  system  would  be  used  to  form  the  "observer  to 
object  vector",  but  the  (X,Y,Z)  coordinates  of  this  point 
are  not  readily  known.  The  observer  location  (X1,Y1,Z1)  is 
not  used  because  it  is  behind  the  image  plane  and  will 
result  in  some  panels  being  viewed  which  are  behind  the 
image  plane.  The  point  (X2,Y2,Z2)  is  used  as  an 
approximation  of  the  origin  of  the  image  plane.  If  the 
angle  between  the  LOS  vector  and  the  object  vector  is 
greater  than  90  degrees,  the  point  lies  behind  the  image 
plane  and  will  not  be  projected  onto  the  image  plane. 
Figure  4-7  illustrates  the  orientation  of  the  vectors  and 
the  image  plane.  The  cosine  of  the  angle  between  the  two 
vectors  is  again  found  by  using  the  dot  product  method. 

If  the  angle  is  less  than  90  degrees,  the  cosine  will 
be  a  positive  number.  If  the  angle  is  greater  than  90 
degrees  the  cosine  will  be  negative;  this  indicates  the 
point  is  not  in  the  hemisphere  in  front  of  the  observer. 
If  the  point  is  not  visible  in  the  forward  hemisphere  then 
the  point  is  flagged  as  hidden.  Generally  each  point  is 
associated  with  up  to  six  triangles;  marking  each  point 
as  hidden  will  prevent  attempting  to  draw  triangles  which 
may  be  partially  behind  the  image  plane. 
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Figure  4-7  Image  Plane  Object  Vector 


The  points  which  are  in  the  forward  hemisphere  are 
projected  onto  the  image  plane  using  the  relation 
developed  in  Chapter  Two,  equation  2.6  which  is  repeated 
here  for  clarity: 


1 
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Expanding  to 

three  separate  equations 

gives: 

*A  *  XB  -  k  [ 

M11 ( XB-XF) 

+  Mi2(YB_YF) 

+ 

M13(2B_2F)  J 

(4.10) 

y»  -  »b  *  k  [ 

H21(XB“XF) 

♦  H22(YB_YF) 

+ 

M23(ZB_ZF)  ] 

(4.11) 
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n* 

B 

*• 

m31 <xB_XF* 

+  M32 ( yb“YF^ 

♦ 

H33  ( 2g-Zp)  ] 

(4.12) 

Then  dividing  (4.10)  and 

(4.11)  by  (4 

.12)  yields: 

r  Hu  (Xg-Xp) 

+  M12(^B~yF) 

+ 

M13(ZB“2F) 

(4.13) 

*  *0  f 

H31(Xg-Xp) 

♦  M32(Vg-Vp) 

+ 

H33  (  2g-Zp) 

M2i(XB-Xp) 

*  m22(YB"'YF) 

+ 

M23  *  ZB“ZF) 

(4.14) 

y  -  y0  -  -£ 

m31(XB“XF) 

+  M32<YB“YF) 

+ 

H33(ZB~ZF) 

Equations  (4.13)  and  (4.14)  are  the  relations  used  to 
transform  the  elevation  data  points  to  the  image  plane. 
They  allow  projecting  all  elevation  points  in  the 
hemisphere  forward  of  the  observer  at  one  time.  Altering 
the  size  of  the  image  plane  will  not  require  the  elevation 
points  to  be  projected  a  second  time.  [Ref.  4,  p.  142] 

D.  AFFINE  IMAGE  PLANE  TO  SCREEN  TRANSFORM 

The  image  plane  coordinates  must  be  transformed  to 
the  screen  coordinate  system.  Figure  4-8  illustrates  the 
image  plane  and  the  screen  coordinate  systems.  The  maximum 
frame  size  values  of  the  x  and  y  axis  in  the  image  plane 
system  are  determined  by  the  desired  magnification  or 
field  of  view  and  the  y-scale  factors. 

In  general,  the  affine  transform  will  map  between  any 
2-D  coordinate  systems  allowing  for  a  rotation  of  the 
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Figure  4-8  Image  to  Screen  Coordinates 


axis,  horizontal  and  vertical  scale  changes,  two 
translations  and  a  non-perpendicularity  of  the  axis.  The  « 

general  form  of  the  affine  transform  is  [Ref.  4,  p.  593]: 


X1 

1 

'  b2 

-bl 

'  x2-  C1 

.  *1  J 

alb2  ~  a2bl 

.  -a2 

al  . 

.  y2~  C2  . 

where:  =  Sx  (cos  P  -£  sin  p) 

a2  =  Sx  (s^n  ^  +  S  cos 
bl  =  sy  (_sin 
b2  =  Sy  (cos  p) 

C  represents  a  non-perpendicularity  of  the  axis 
P  represents  a  rotation  of  the  axis 
=  x  translation  of  origin 
C2  =  y  translation  of  origin 
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Starting  with  a  magnification  factor  of  1.0  and  a  y- 
scale  factor  of  1.2,  the  x  and  y  frame  sizes  are  35000 
and  29000  respectively.  These  values  must  be  mapped  onto  a 
screen  coordinate  system  which  ranges  from  1-512  in  both 
the  x  and  y  direction.  There  are  two  scale  changes,  two 
translations,  no  rotations  and  no  non-perpendicularities 
involved  with  this  transformation.  This  results  in  sigma 
equal  to  zero  and  beta  equal  to  zero.  Substituting  these 
values  results  in: 


(4.17) 


The  scale  factors  are  formed  at  runtime  to  allow  the 
changes  in  magnification  (field  of  view)  and  elevation 
scale.  Figure  4-9  is  psuedocode  which  implements  the  2-D 
transform.  Note  that  the  variables  Al,  A2 ,  Bl,  B2  are  in 
terms  of  the  image  plane  frame  size.  This  approach  allows 
the  user  to  change  the  magnification  or  elevation  scaling 
interactively  at  this  point.  Also  note  that  each  elevation 
point  is  checked  for  non-visible  status  prior  to  the 
transform.  This  checking  saves  execution  time  by  avoiding 
points  not  in  the  field  of  view. 
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BEGIN 

Cl  =  image  x  axis  maximum 

C2  =  image  y  axis  maximum 

A1  =  image  x  axis  frame  size  /  screen  x  axis  maximum 

B2  =  image  y  axis  frame  size  /  screen  y  axis  maximum 

DO  IR  *  1  to  Number  of  elevation  pts. 

IF  (point  .  EQ.  visible)  THEN 
xscreen  =  (ximag  -  Cl)  /  A1 

yscreen  =  (yimag  -  C2)  /  B2 

END  IF 

END  DO 


END 


Figure  4-9  Affine  Transform  Psuedocode 


E.  HIDDEN  SURFACE  REMOVAL 

In  order  to  properly  display  the  perspective  view,  the 
surfaces  hidden  behind  front  surfaces  must  be  dealt  with. 
Earlier  work  utilized  a  Z-buffer  to  perform  the  hidden 
surface  removal  [Ref.  2,  p.  41].  The  Z-buffer  method 
utilizes  a  pixel  by  pixel  depth  buffer.  Each  pixel 
position  on  the  screen  corresponds  to  a  location  in  the 
buffer.  The  buffer  is  initialized  to  a  maximum  depth,  then 
prior  to  drawing  each  pixel  the  depth  is  compared  with  the 
depth  already  stored  in  the  buffer.  If  the  pixel  depth  is 
less  than  the  depth  stored  in  the  buffer,  then  the  pixel 
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is  drawn  and  the  new  depth  is  stored  in  the  depth  buffer. 
Earlier  work  by  Coleman  pointed  out  the  poor  performance 
of  this  technique  when  executed  in  software  [Ref.  2,  p. 
55].  Several  other  methods  of  hidden  surface  removal  were 
investigated;  each  will  be  discussed  briefly  and  the 
disadvantages  listed. 

First,  the  method  of  ’’bounding  rectangles”  was 
investigated.  In  this  method  a  rectangle  is  formed  around 
each  triangle.  This  rectangle  is  as  small  as  possible  and 
is  constructed  from  the  three  vertices  forming  the 
triangle.  The  (x,y)  coordinates  of  each  vertex  are  checked 
and  the  minimum  and  maximums  of  both  the  x  and  y 
coordinates  form  the  boundaries  of  the  rectangle.  Now  the 
vertices  of  all  other  triangles  are  compared  to  the 
rectangle.  If  the  (x,y)  coordinates  of  a  vertex  lie  within 
the  boundary  of  the  rectangle,  the  triangle  is  marked  as 
overlapping  and  the  Z-buffer  routine  is  used  to  handle  the 
hidden  surface  removal.  If  the  triangle  does  not  overlap 
any  other  triangles,  it  is  drawn  to  the  screen.  Originally 
the  intent  was  to  reduce  the  number  of  panels  which  needed 
to  be  drawn  using  the  Z-buffer  method.  This  was  not 
achieved  because  the  panels  are  triangular  vice 
rectangular.  Several  experiments  showed  that  no  triangles 
passed  the  bounding  rectangle  test  and  consequently  no 
savings  in  execution  time  was  realized.  [Ref.  9,  p.  246] 
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A  second  method  utilizes  a  test  function  for  each  side 
of  a  triangle.  A  test  function  is  formed  for  each  side  of 
a  triangle.  This  test  function  is  formed  from  the  equation 
of  the  line  describing  the  side.  The  vertices  of  all  other 
triangles  are  checked  against  the  test  function.  This 
results  in  nine  comparisons  for  every  triangle  checked. 
The  calculation  time  is  excessive,  also  the  results  were 
similar  to  the  results  of  the  bounding  box  test.  Very  few 
triangles  passed  the  test  and  the  Z-buffer  method  had  to 
be  used  on  the  majority  of  panels,  with  no  savings  in 
execution  time. [Ref.  9,  p.  24] 

Another  approach,  known  as  the  "Painter's  algorithm", 
is  not  hidden  surface  removal  at  all.  [Ref.  8,  p.  265] 
Here  the  panels  are  drawn  to  the  screen  in  sequence  from 
background  to  foreground.  The  panels  which  are  hidden  from 
view  are  covered  with  the  panels  which  are  closer  to  the 
observer.  This  is  the  approach  used  in  this  software.  Each 
panel  is  sorted  into  order  based  on  the  maximum  distance 
of  the  vertices  from  the  observer  location.  The  choice  of 
sorting  algorithms  has  a  large  affect  on  the  efficiency  of 
this  method.  Originally,  the  routine  incorporated  a  simple 
bubble  sort.  The  cost  in  performance  is  not  noticeable 
with  a  small  number  of  panels,  but  as  the  resolution  of 
the  bathymetric  data  increases,  the  number  of  panels 
increases  as  shown  in  equation  3.1  which  is  repeated  here. 
#  panels  =  ((#rows  -  1) *2 )  *  (#cols  -  1)  (4.18) 
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The  poor  performance  of  the  bubble  sort  becomes  apparent 
as  the  number  of  panels  increases  beyond  500.  Two 
alternate  sorting  routines  were  used  with  this  routine, 
the  shell  sort  and  the  quick  sort.  Table  4-1  is  a  table 
which  shows  the  measured  performance  of  these  sort 
routines.  Clearly  for  large  numbers  of  panels  the  quick 
sort  is  the  most  efficient.  Based  on  these  results  the 
quick  sort  was  selected  for  use  in  the  hidden  surface 
routine.  Figure  4-11  is  psuedocode  which  describes  the 
hidden  surface  removal  procedure. 


TABLE 

4-1  SORT 

Time 

PERFORMANCE 

in  seconds 

#  panels 

Bubbl e 

Shell 

Quick 

500 

2 

<1 

<1 

1000 

8 

<1 

<1 

2000 

33 

<1 

<1 

3000 

80 

1 

1 

5000 

210 

2 

1 

10000 

900 

5 

1.5 

20000 

— 

14 

3 

50000 

— 

56 

9 
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BEGIN 


(*  create  sort  key  by  counting  the  panels  in  front  of 
the  image  plane.  Depth_key(L)  contains  the  panel  ID 
number  *) 

DO  IR  -  1  to  NPLANES 

IF  (  panel  .NOT.  behind  field  of  view)  THEN 
INCREMENT  COUNT 
Depth_key( COUNT)  =  IR 
END  IF 

END  DO 

(*  calc  maximum  depth  of  each  visible  plane  *) 

DO  L  =  1  to  COUNT 
IR  =  Depth_key(L) 

PL_AR ( IR, 5) = 

MAX ( DEPTH (node_a) , DEPTH (node_b) , DEPTH (node_c) ) 

END  DO 

(*  call  sort  routine  to  sort  the  panels  on  the 
maximum  depth*) 

CALL  QUICKSORT 

END 


Figure  4-11  Psuedocode  for  Hidden  Surface  Removal 
by  Sorting 


Upon  completion  of  the  hidden  surface  removal  routine, 
the  screen  is  erased  and  the  panels  are  plotted.  The 
vertices  for  each  panel  along  with  the  shading  are  stored 
in  array  PL_AR(*,*).  The  shading  value  is  clipped  to  a 
range  of  1-250  vice  the  original  1-256  due  to  the  DEC  VMS 


window  manager.  This  maintains  the  background  color  of  the 
screen  and  the  window  attributes.  The  panels  are  read  in 
sequence  using  the  sort  index.  The  (x,y)  screen 
coordinates  of  each  vertex  and  the  shade  are  passed  to  the 
plot  subroutine  and  the  panel  is  drawn  on  the  screen. 


F.  RESULTS 

Data  cropped  out  of  the  primary  bathymetry  data  file 
covers  an  area  bounded  by  37:00:00  N.  latitude,  126:00:00 
W.  longitude  to  36:00:00  N.  latitude,  125:00:00  W. 
longitude.  Figure  4-12  is  the  digitized  image  of  the 
mosaic  used  for  the  image  data.  The  first  synthesized  view 
Figure  4-13,  is  from  an  observer  location  of  36:0:0  N. 
125:30:0  W. ,  4000  meters  below  sea  level  and  heading 
equal  to  000  degrees.  The  y-scale  exaggeration  is  six  and 
the  magnification  factor  is  one.  Figures  4-14  a  to  c  are 
displays  of  the  same  area  from  different  observer 
locations.  Figures  4-l4a  and  4 -14b  are  offset  from  the 
original  location  east  and  west  by  25  minutes  of 
longitude.  Figure  4-14c  is  from  the  same  location  as  the 
original  but  the  height  is  2000  meters  below  sea  level. 
These  images  display  the  ability  to  "see”  the  ocean  bottom 
in  a  new  fashion.  The  resolution  of  the  input  elevation 
data  was  originally  5  minutes  by  5  minutes.  This  data  was 
interpolated  in  both  latitude  and  longitude  to  a 
resolution  of  1  minute  by  1  minute  in  order  to  achieve  a 
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better  shading  rendition.  This  reduces  the  number  of 
pixels  averaged  for  each  triangle  and  improves  the  shading 
of  the  generated  image.  Figure  4-15  is  an  image  using  the 
original  5  minute  by  5  minute  resolution  data.  Notice 
that  the  size  of  the  panels  is  larger  and  the  shading 
rendition  is  much  more  coarse.  In  addition  to  affecting 
the  shading,  the  resolution  of  the  elevation  data  directly 
affects  the  physical  shaping  of  objects  in  the  terrain. 
This  program  does  not  utilize  any  method  of  curve  fitting 
between  adjacent  elevation  points  and  simply  draws 
straight  lines  between  the  points.  This  results  in 
objects  being  coarsely  approximated  in  the  synthesized 
view. 


Figure  4-12  Mosaic  Image 


Figure  4-13  Obs  Loc:  36:00:00  N  Depth  4000  m 

125:30:00  W  Heading  000  deg 
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Figure  4-l4b  Obs  Loc  36:00:00  N  Depth  4000  m 

125:55:00  W  Heading  020  deg 
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V.  OPERATIONAL  ASPECTS  OF  THE  PROGRAM 


The  software  utilizes  an  input  file  to  load  the  names 
and  geographic  locations  of  the  data  files.  This  input 
file  may  be  created  by  an  ASCII  editor  or  created  within 
the  program  by  following  the  prompts.  The  format  of  the 
input  file  is  shown  in  Figure  5-1. 

If  the  input  file  exists,  the  user  simply  enters  the 
name  when  prompted.  The  program  will  seem  to  pause  at  this 
point  while  reading  in  the  data  files,  forming  the  panels 
and  averaging  the  image  pixels.  When  complete  with  these 
steps,  the  software  will  prompt  for  the  observer  location. 
This  location  is  entered  in  terms  of  observer  latitude  and 
longitude  in  degrees,  minutes  and  seconds.  Southern 
latitudes  and  western  longitudes  are  entered  as  negative 
values.  Also  entered  is  the  elevation  of  the  observer  and 
the  heading  of  the  observer.  Elevation  is  in  meters  with 
respect  to  sea  level  with  values  greater  than  sea  level 
entered  as  positive  values  and  elevations  below  sea  level 
entered  as  negative  values.  The  heading  is  entered  as  a 
value  from  0-360  degrees  where  0  deg.  is  north,  90  deg.  is 
east,  180  deg.  is  south,  270  deg.  is  west  and  360  deg.  is 
north  again. 
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Elevation  data  file  name 


#  of  rows  of  elevation  data 

#  of  columns  of  elevation  data 
Image  data  file  name 

#  of  rows  of  image  data 

#  of  columns  of  image  data 

Latitude/Longitude  of  Northwest  corner  of  image 
data  in  (deg, min, sec) 

Latitude/ Longitude  of  Southeast  corner  of  image 
data  in  (deg, min, sec) 

Title  for  the  screen  image 

Figure  5-1  Input  File  Format 


At  this  point  the  software  will  draw  the  perspective 
view  on  the  screen.  It  takes  40  seconds  to  produce  one 
image  using  61  x  61  elevation  points  and  512  x  512  image 
points.  After  the  image  is  drawn  a  menu  is  presented.  The 
operator  may  alter  the  magnification  or  field  of  view, 
change  the  elevation  scale  exaggeration,  enter  a  new 
observer  location,  save  the  image  or  exit  the  program. 

A.  ALTER  MAGNIFICATION  OR  FIELD  OF  VIEW 

The  perspective  view  transformation  simulates  viewing 
through  a  camera  viewfinder.  The  observer  selects  the 
location  and  points  the  camera,  and  the  image  is  formed  in 
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the  viewfinder.  The  magnification  or  field  of  view  is 
changed  by  changing  the  focal  length  of  the  lens. 

The  magnification  factor  of  the  initial  image  is  one. 
The  values  of  the  image  plane  x-axis  and  y-axis  maximums 
are  initially  set  to  35000  and  29000  respectively.  This 
loosely  corresponds  to  a  35mm  x  29mm  frame  size  of  a  35mm 
camera.  The  focal  length  corresponds  to  an  initial  focal 
length  of  50mm.  The  field  of  view  is  given  by: 

,  f  xframe  ^ 

Field  of  View  -  2  tan"1 £ocua)  J  (5.1) 

Using  the  values  above  results  in  a  field  of  view  of  36.8 
degrees.  This  field  of  view  is  approximately  equal  to  the 
field  of  view  from  a  35mm  camera  with  a  50mm  "normal" 
lens.  This  field  of  view  corresponds  to  a  magnification 
factor  of  one.  Table  5-1  shows  the  field  of  view  for 
various  magnifications. 

When  changing  the  magnification  on  a  camera,  the 
standard  method  is  to  change  the  focal  length  of  the  lens. 
This  method  presents  a  disadvantage  here  in  that  it 
requires  the  complete  3-D  to  3-D  transform  to  be 
performed  for  every  magnification  change.  By  changing  the 
frame  size  instead  of  the  focal  length,  the  new  image  is 
generated  faster  because  the  change  is  made  during  the 
affine  image  plane  to  screen  plane  transformation.  The 
magnification  changes  are  entered  as  values  with  respect 
to  the  normal  magnification. 


TABLE  5-1  MAGNIFICATION  VERSUS  FIELD  OF  VIEW 

Magnification 

Field  of  View  (deg) 

.5 

69.9 

1.0 

38.5 

1.5 

26.2 

2.0 

19.8 

2.5 

15.9 

3.0 

13.3 

3.5 

11.4 

4.0 

10.0 

4.5 

8.9 

5.0 

8.0 

B.  CHANGING  THE  ELEVATION  SCALING 

Altering  the  y-axis  scale  allows  exaggeration  or 
reduction  of  the  y-scale.  Views  of  the  terrain  with  the 
elevation  scaling  set  at  1:1  result  in  views  which  have 
very  little  appearance  of  height.  By  modeling  the  frame 
size  as  a  35mm  camera  frame,  a  1.2:1  scaling  is 
introduced.  Research  conducted  by  the  U.S.  Army  Research 
Institute  for  the  Behavioral  and  Social  Sciences  indicated 
that  vertical  scaling  ranging  from  1.25:1  to  1.50:1 
presented  the  most  natural  views.  [Ref.  3,  p.  37] 

The  elevation  scaling  is  changed  in  a  fashion  similar 
to  the  magnification  as  a  value  with  respect  to  the  1:1 
scale. 
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C.  ENTER  A  NEW  OBSERVER  LOCATION 

This  option  allows  the  operator  to  enter  a  new 
observer  location.  The  current  settings  for  magnification 
and  elevation  scaling  are  maintained. 

D.  SAVING  AN  IMAGE 

Selecting  this  option  allows  multiple  images  to  be 
saved  to  a  disk  file.  When  the  option  is  selected  the 
first  time,  it  will  prompt  for  the  file  name  for  the 
images  to  be  stored  in.  Two  data  files  are  created,  *.IMG 
and  *.SIZ.  If  the  filename  extension  is  specified  then  it 
is  used  instead  of  the  .IMG  extension.  Images  will  be 
saved  to  *.IMG  file  for  the  entire  session  each  time  the 
(S)ave  option  is  selected.  The  entire  set  of  images  may  be 
viewed  sequentially  using  a  program  named  PLAYBACK. FOR. 
This  will  display  all  the  images  stored  in  the  image  file 
at  a  rav.e  of  one  image  per  second.  Another  program, 
UIST0512 . FOR  will  convert  a  single  frame  from  the  saved 
images  to  a  512  x  512  row-column  format,  with  512  fixed 
length  records  and  512  bytes  per  record.  The  *.SIZ  file 
contains  the  buffer  size  data  required  by  the  playback 
routine. 
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VI.  CONCLUSIONS 


A.  GENERAL 

The  goal  of  this  thesis  was  creating  the  3-D 
perspective  view  of  ocean  bottom  terrain  and  applying 
shading  or  color  information  from  an  image  file.  This  goal 
was  achieved  and  the  ability  to  Mseew  the  terrain  in  this 
fashion  holds  promise  for  future  work  in  many  fields. 

The  program  presented  here  combines  terrain  elevation 
data  with  shading  information  from  a  second  data  source. 
This  fusion  of  bathymetry  and  imagery  data  results  in  a 
3-D  perspective  view  of  the  ocean  bottom  in  a  fashion 
which  appears  natural  to  the  observer.  This  type  of 
display  system,  which  combines  multiple  sources  of  data 
into  a  form  that  is  easier  to  comprehend,  has  potential 
uses  in  mapping,  geologic  exploration  and  weapon  system 
displays.  It  has  direct  use  in  any  geographic  information 
system  where  the  data  stored  consists  of  (x,y)  coordinates 
and  an  attribute.  This  program  could  merge  any  two  layers 
of  information  from  such  a  system  and  display  a  view  which 
is  easier  to  interpret.  The  data  utilized  in  this  project 
was  not  a  part  of  such  a  system,  but  easily  could  have 
been.  In  this  case  the  first  layer  of  data  would  be  the 
(x,y)  location  and  the  elevation  with  respect  to  sea  level 
and  the  second  layer  of  data  would  be  the  (x,y)  location 
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and  the  shading  value.  Combining  the  layers  of  data  from 
multiple  sources  is  a  form  of  sensor  fusion. 


B.  LIMITATIONS 

There  are  several  limitations  in  this  program.  First, 
the  image  plane  orientation  is  fixed  in  the  vertical 
direction  with  respect  to  the  observer  heading  and 
location.  This  implies  that  the  perspective  view  formed 
simulates  the  view  from  a  vehicle  in  level  flight. 
Allowing  the  image  plane  to  deviate  from  vertical  would 
allow  the  observer  to  look  down  on  the  terrain  and 
generate  views  not  presently  permitted.  Second,  the 
shading  values  are  clipped  to  250  vice  255  shades.  This  is 
done  in  order  to  maintain  five  separate  entries  in  the 
virtual  color  map  which  correspond  to  the  system  colors.  A 
third  problem  exists  with  the  geographic  areas  of  the 
input  data.  The  conversion  from  latitude/longitude 
coordinates  to  geo-rectangular  coordinates  is  not 
implemented  at  the  hemisphere  boundaries  such  as  at  the 
equator  (0  deg.  N.),  or  at  0  deg.  longitude.  The  input 
data  must  be  completely  above  or  below  the  equator  and 
similarly  completely  east  or  west  of  0  deg.  longitude. 


C.  PERFORMANCE 

Forming  the  3-D  perspective  views  requires  the  input 
data  to  be  transformed  from  one  3-D  coordinate  system  to 
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the  image  plane  system  and  ultimately  to  the  screen 
coordinate  system.  As  discussed  in  Chapter  Three,  the  size 
of  the  input  data  files,  specifically  the  elevation  data 
will  affect  the  performance  of  the  program.  Table  6-1  is 
a  breakdown  of  the  amount  of  time  required  to  complete 
each  of  the  major  sections  of  the  program.  The  overall 
performance  is  satisfactory  and  breaks  down  to  40  seconds 
to  complete  a  view  after  the  observer  location  is  entered. 
The  plotting  of  the  view  on  the  screen  is  89%  of  this 
time,  while  the  calculation  involved  accounts  for  only  11% 
of  the  time.  Clearly  the  plotting  time  must  be  improved 
prior  to  spending  an  appreciable  amount  of  time  in  the 
calculation  subroutines.  The  length  of  time  required  by 
the  plotting  routine  is  related  to  the  number  of  triangles 
being  passed  to  the  graphics  hardware.  The  communication 
involved  in  passing  the  vertices  of  each  panel  to  the 
graphics  processor  is  the  bottleneck.  Improvements  may  be 
achieved  in  this  area  by  writing  a  driver  to  directly  load 
the  appropriate  values  to  the  graphics  hardware  and 
bypassing  the  overhead  of  the  software  routines  presently 
utilized. 

The  combination  of  the  perspective  transformation  and 
the  hidden  surface  removal  routines  account  for  95%  of  the 
11%  calculation  time  as  seen  in  Table  6-1.  Other  methods 
of  achieving  slight  performance  increases  may  include  a 
hardware  implementation  of  the  perspective  transformation 


55 


TABLE  6-1  COMPUTATION  TIME 
Subsection  Resolution 


61  x  61  13  x  13 

(time  in  seconds) 

Set  up  graphics  .80  .53 

User  input  -  - 

Read  elev.  file  .25  .11 

Read  image  file  10.22  10.69 

Form  triangles  3.3  .3 

Average  gray  shade  6.87  2.54 

Observer  location  -  - 

Form  [M]  matrix  .01  .01 

[M]  multiply  2.51  .12 

Affine  transform  .23  .01 

Hidden  surface  removal  1.79  .07 

Plotting  35.24  1.76 

Total  elapsed  after  -  - 

Observer  location  39.78  1.97 

Calculation  time 

percentage  11.4%  10.6% 


matrix  multiplies  and  a  more  efficient  hidden  surface 
removal  routine.  Altering  the  hidden  surface  removal 
routine  most  likely  will  result  in  more  calculation  time 
required  in  that  area,  but  may  result  in  substantial 
savings  in  the  plotting  time.  The  approach  for  hidden 
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surface  removal  used  here  does  not  actually  remove  hidden 
triangles  but  draws  them  in  order  from  the  background  to 
foreground.  By  altering  the  hidden  surface  routine  so  that 
hidden  panels  are  not  plotted  at  all,  a  savings  may  be 
made  in  the  time  required  to  plot  the  view. 

D.  FUTURE  WORK 

This  program  is  a  first  step  in  the  direction  of 
combining  various  layers  of  data  into  a  single  form.  There 
are  several  areas  which  are  being  considered  for  further 
work. 

The  type  of  data  merged  onto  the  elevation  grid  is  not 
limited  to  the  imagery  data  used  here.  Magnetic  anomaly 
data  is  available  and  could  also  be  merged  onto  the  grid. 
This  would  permit  the  observer  to  make  correlations 
between  the  terrain  and  the  magnetic  response. 

The  original  effort  in  this  area  began  with  aerial 
photographs  and  surface  terrain.  Applying  the  terrain 
elevation  data  and  an  aerial  photograph  to  this  program 
would  broaden  the  versatility  of  the  program  greatly. 

Combining  two  layers  of  data  is  a  starting  point,  the 
ability  to  combine  three  or  four  layers  onto  the 
elevation  grid  using  colors  to  distinguish  the  layers  may 
be  possible. 

In  order  to  make  this  program  truly  useful,  a  user 
interface  must  be  designed.  As  a  minimum  this  requires  a 


method  of  entering  observer  locations  based  on  the  current 
location  and  perspective  view.  A  particularly  effective 
user  interface  would  involve  using  a  pointer  to  select  an 
area  on  the  screen  and  tell  the  program  "I  wish  to  go 
there".  This  would  permit  "driving  around"  an  area  of 
interest . 

Finally,  a  method  of  loading  data  from  the  storage 
device  and  merging  it  with  the  currently  loaded  data  is 
required.  This  would  remove  the  boundaries  from  the  image 
displayed  on  the  screen  so  that  when  the  observer  location 
is  near  the  edge  of  the  current  area,  the  adjacent  areas 
will  be  loaded  and  displayed. 


APPENDIX  A. 


This  appendix  lists  the  different  routines  of  the 
program  by  subroutine  name  and  briefly  lists  the  inputs 
and  outputs  of  the  subroutine.  Also  given  are  the  names  of 
the  subroutines  which  interact  with  each  subroutine. 

1.  S0NAR3D 

A.  Function 

Main  routine,  calls  other  supporting  routines  to 
complete  the  3-D  perspective  view 

B.  Input 
None 

C .  Output 
None 

D.  Calling  routine 
None 

E.  Called  routine 
INPUT 
READ_ELEV 

RE A D_ I MAGE 
TER_XYZ 
IM_REFAVG 
OBS_LOC 
M  ORIEN 
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NEW  IJ 


XY2IJ 

HIDDENJSURF 

CLRSCRN 

2.  Subroutine  INPUT 

A.  Function 

Reads  in  or  creates  the  initial  data  required  for 
the  elevation  data  file  and  the  image  data  file. 

B.  Input 

NAME  -  name  of  input  data  file 

ELFILE  -  elevation  data  file  name 

EL_SIZ  -  elevation  data  number  of  columns 

EL_REC  -  elevation  data  number  of  rows 

IMFILE  -  image  data  file  name 

IM_SIZ  -  image  data  number  of  columns 

IM_REC  -  image  data  number  of  rows 


I  LAD, 

I LAM, ILAS  -  northwest 

latitude 

of 

image 

file 

ILOD, 

ILOM, ILOS  -  northwest 

longitude 

of 

image 

file 

FLAD, 

FLAM,FLAS  -  southeast 

latitude 

of 

image 

file 

FLAD, 

FLAM, FLAS  -  southeast 

longitude 

of 

image 

file 

TITLE 

-  50  character  title  placed 

on 

the 

image 

screen 

C .  Output 

Same  as  input 

D.  Calling  routines 
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S0NAR3D 


E.  Called  routines 
CLRSCRN 

3 .  Subroutine  READ_ELEV 

A.  Function 

Read  in  the  elevation  data  file  from  disk 

B.  Input 
None 

C .  Output 

ILATD, ILATM, ILATS  -  southwest  latitude  of  the 

elevation  data  file 

ILOND, ILONM, ILONS  -  southwest  longitude  of  the 

elevation  data  file 

D_LAT  -  latitude  grid  size  in  seconds 

D_LONG  -  longitude  grid  size  in  seconds 

IENDM  -  number  of  rows  of  the  elevation  data 

IENDN  -  number  of  columns  of  the  elevation  data 

RELEV (*,*)  -  array  containing  the  elevation  data 

points 

D.  Calling  Routines 
S0NAR3D 

E.  Called  routines 
None 
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4 .  Subroutine  READ_IMAGE 

A.  Function 

Read  in  the  image  data  file  from  disk 

B.  Input 

IM_SIZ  -  number  of  columns  of  the  image  data 
IM_REC  -  number  of  rows  of  the  image  data 

C .  Output 

IMAGE (*,*)  -  array  containing  the  image  data 

D.  Calling  routines 
S0NAR3D 

E.  Called  routines 
None 

5.  Subroutine  TER_XYZ 

A.  Function 

Converts  each  elevation  point  to  it's  geo- 
rectangular  XYZ  coordinates  and  calculates  the 
corresponding  image  row/column  coordinates  for  each 
point 

B.  Input 

LATD, LATM, LATS  -  reference  elevation  latitude 
ILOND, LONM, LONS  -  reference  elevation  longitude 
D_LAT  -  latitude  grid  size  in  seconds 
D_LONG  -  longitude  grid  size  in  seconds 
RELEV ( * , * )  -  array  holding  the  elevation  data 
IENDM  -  number  of  rows  of  the  elevation  data 
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IENDN  -  number  of  columns  of  the  elevation  data 
ILAD, ILAM, ILAS  -  reference  northwest  image  latitude 
ILOD , ILOM, ILOS  -  reference  northwest  image 
longitude 

FLAD, FLAM, FLAS  -  reference  southeast  image  latitude 
FLOD ,  FLOM ,  FLOS  -  reference  southeast  image 
longitude 

IM_SIZ  -  number  of  columns  of  the  image  data 
IM_REC  -  number  of  rows  of  the  image  data 

C .  Output 

IA(*)  -  array  containing  the  image  column 

coordinates 

JA ( * )  -  array  containing  the  image  row  coordinates 
XYZ ( * , 3 )  -  array  holding  the  (X,Y,Z)  coordinates  of 
each  elevation  point 

D.  Calling  routines 
SQNAR3D 

E.  Called  routines 
C0NV2SEC 
DMS2XYZ 

6.  Subroutine  IM_REFAVG 

A.  Function 

This  routine  forms  the  panel  array  containing  the 
nodes  of  the  triangles  and  the  average  gray  shade. 

B.  Inputs 
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IA ( * )  -  array  containing  the  image  column 

coordinates 

JA(*)  -  array  containing  the  image  row  coordinates 
IENDM  -  numbe  r  of  rows  of  the  elevation  data 
IENDN  -  number  of  columns  of  the  elevation  data 
IMAGE (*,*)  -  array  containing  the  image  data 

C .  Outputs 

PL_AR(*,5)  -  array  containing  the  nodes  of  the 

triangle  and  the  gray  shade 

D.  Calling  routines 
S0NAR3D 

E.  Called  routines 
None 

7 .  Subroutine  OBSLOC 

A.  Function 

Calculates  the  observer  location  in  (X,Y,Z) 
coordinates  given  the  latitude  and  longitude. 

B.  Input 

LATD , LATM , LATS  -  Latitude  of  observer  location 
LON D, LONM, LONS  -  Longitude  of  observer  location 
HEIGHT  -  Observer  elevation  in  meters 
CTS  -  Observer  course  in  degrees 

C.  Outputs 

OBS_LOC ( * )  -  array  holding  observer  location 

parameters 
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X1,Y1,Z1  -  Geo-rectangular  coordinates  of  observer 
location 


X2,Y2,Z2  -  Geo-rectangular  coordinates  of  second 

observer  point  along  LOS 

D.  Calling  routines 
S0NAR3D 

E.  Called  routines 
DMS2XYZ 

8.  Subroutine  M_ORIEN 

A.  Function 

Determine  the  orientation  of  the  image  plane  and 
calculate  the  [M]  matrix  parameters 

B.  Inputs 

X1,Y1,Z1  -  Geo-rectangular  coordinates  of  observer 
location 

X2,Y2,Z2  -  Geo-rectangular  coordinates  of  second 

observer  point  along  LOS 

C.  Outputs 

M(3,3)  -  [M]  matrix  parameters 

D.  Calling  routines 
S0NAR3D 

E.  Called  routines 


None 
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9 .  Subroutine  NEW_IJ 

A.  Function 

Calculate  the  image  plane  coordinates  of  all  object 
points  which  are  in  front  of  the  field  of  view. 
Also  mark  as  HIDDEN  any  nodes  behind  the  image 
plane 

B.  Inputs 

X1,Y1,Z1  -  Geo-rectangular  coordinates  of  the 

observer  location 

X2,Y2,Z2  -  Geo-rectangular  coordinates  of  the 

second  observer  point  along  LOS 

ITOT  -  Total  number  of  elevation  points 
* 

XYZ ( * , 3 )  -  array  holding  the  (X,Y,Z)  coordinates  of 
each  elevation  point 
focus  -  Focal  length 
M ( 3 , 3 )  -  [M]  matrix  parameters 

C .  Outputs 

IMAX ( * )  -  x  image  coordinate 
IMAY(*)  -  y  image  coordinate 

DEPTH (*)  -  Distance  from  observer  of  each  elevation 
point 

HIDDEN (*)  -  Flag  for  points  hidden  behind  image 

plane 

D.  Calling  routines 
S0NAR3D 
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E.  Called  routines 
None 

10.  Subroutine  XY2IJ 

A.  Function 

Convert  the  image  (x,y)  coordinates  to  screen  (i,j) 
coordinates 

B.  Inputs 

IMAX ( * )  -  x  image  coordinate 

IMAY ( *)  -  y  image  coordinate 

ITOT  -  Total  number  of  elevation  points 

XIMAMAX  -  image  frame  size  (x  direction) 

YIMAMAX  -  image  frame  size  (y  direction) 

HIDDEN (*)  -  Flag  for  points  hidden  behind  image 

plane 

C .  Outputs 

IA(*)  -  screen  x  coordinate 
JA(*)  -  screen  y  coordinate 

D.  Calling  routines 
SONAR3D 

E.  Called  routines 
AFFIN 

11.  Subroutine  AFFIN 
A.  Function 
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Calculates  the  coefficients  for  the  AFFIN  transform 
from  image  coordinates  to  screen  coordinates 

B.  Inputs 

XIMA_MAX  -  image  frame  size  (x  direction) 

YIMA_MAX  -  image  frame  size  (y  direction) 

C.  Outputs 

A1 ,  A2  ,  B1 ,  B2  ,  Cl  ,C2  -  parameters  for  the  AFFIN 

transform 

D.  Calling  routines 
XY2IJ 

E.  Called  routines 
None 

12.  Subroutine  HIDDENJSURF 

A.  Function 

Remove  hidden  surfaces  and  plot  to  screen 

B.  Inputs 

IA(*)  -  screen  x  coordinate 
JA(*)  -  screen  y  coordinate 

IENDM  -  Number  of  rows  of  the  elevation  data 
IENDN  -  Number  of  columns  of  the  elevation  data 
DEPTH  (*)  -  Distance  from  the  observer  of  each 

elevation  point 

PL_AR(*,5)  -  array  containing  the  nodes  of  the 

triangle  and  the  gray  shade 
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HIDDEN (*)  -  Flag  for  points  hidden  behind  the  image 
plane 

VD_ID  -  Virtual  display  identifier 
WD_ID  -  Window  identifier 

C .  Output 
None 

D.  Calling  routine 
SONAR3D 

E.  Called  routine 
QUICKSORT 
UISDC$ERASE 
UIS$SET_WRITING_INDEX 
UISDC$PLOT 

13.  Subroutine  QUICKSORT 

A.  Function 

Sort  the  panels  based  on  maximum  distance  from  the 
observer 

B.  Input 

ARRAY(*,5)  -  array  to  sort 
KEY (*)  -  index  to  main  array 
COUNT  -  number  of  elements  to  sort 

C .  Output 

ARRAY(*,5)  -  original  array 

KEY ( * )  -  index  to  main  array  (revised  order) 
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D.  Calling  routine 
H I DDEN_SURF 

E.  Called  routine 
None 

14 .  Subroutine  CLRSCRN 

A.  Function 

Clear  the  screen 

B.  Input 
None 

C .  Output 
None 

D.  Calling  routine 
SONAR3D 

INPUT 

E.  Called  routines 
None 

15.  Subroutine  C0NV2SEC 

A.  Function 

Convert  DEG, MIN, SEC  coordinates  to  seconds 

B.  Input 

DEG, MIN, SEC  -  Latitude  or  longitude  in  deg, min, sec 
format 

C .  Output 

SEC  -  Total  number  of  seconds 


I 
1 


70 


D.  Calling  routine 
TER_XYZ 

E.  Called  routine 
None 

1  .  Subroutine  DMS2XYZ 

A.  Function 

Convert  DMS  data  to  (X,Y,Z)  coordinates 

B.  Input 

LATD, LATM, LATS  -  Latitude  in  deg, min, sec  format 
LOND , LONM , LONS  -  Longitude  in  deg, min, sec  format 
HEIGHT  -  elevation  with  respect  to  sea  level  in 
meters 

C .  Output 

X, Y, Z  -  (X, Y, Z)  coordinates  of  input  point 

D.  Calling  routine 
TER_XYZ 
OBS_LOC 

E.  Called  routine 
None 
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APPENDIX  B.  PROGRAM  LISTING 


PROGRAM  S0NAR3D 

************************************************************ 

THIS  PROGRAM  CONSTRUCTS  THE  ELEVATION  FILE  AND  IMAGE 
FILE  THAT  IS  REQUIRED  BY  THE  3-D  PROGRAM.  MAXIMUM 
ELEVATION  ARRAY  SIZE  IS  70  ROWS  X  70  COLUMNS.  MAXIMUM 
IMAGE  SIZE  IS  1024  X  1024. 

************************************************************ 
IMPLICIT  INTEGER  (A-Z) 

get  graphics  libraries 

INCLUDE  1  SYS$LEBRARY :  UISUSRDEF 1 
INCLUDE  '  SYS$LLBRARY :  UISENIRY ' 

INCLUDE  'SYS$LLBRAKY:UISMSG' 

C 

C  the  following  variables  define  the  elevation  file  maximum  size 
C  and  the  image  file  maximum  size.  In  order  to  simplify  altering 
C  of  these  values,  they  have  been  grouped  together.  Change  the 
C  values  in  the  DATA  statement,  and  in  the  definition  statements 
C  which  immediately  follow. 

C 

OCMM3N  /ELEV/RCWJELEV,  OOL_ELEV,TOT_ELEV,NFIANES 
COMMON  /IMAGE/RCW_IMAGE ,  COL_IMAGE 

DATA  RCW_ELEV,  OOL_ELEV,  TOT_ELEV ,  NPLANES/70 ,70,4900, 10000/ 
DATA  ROW_IMAGE,OOL_IMAGE/1024 , 1024/ 

C 

BYTE  IMAGE (1024, 1024) 

INTEGER  HIDDEN(4900)  ,IA(4900)  ,JA(4900) 

INTEGER  FL_AR  ( 10000 , 5 ) ,  DEFTH_KEY  ( 10000 ) 

REAL  REIEV(70,70) 

REAL*8  DEPIH(4900) ,XYZ(4900,3) ,IMAX(4900) ,IMAY(4900) 

C 

CHARACTER  ELFIIE*  13 ,  IMFII£*13 ,  FII£_NAME*13 ,  ANS*1 ,  TTTLE*50 
C 

INTEGER  IENEN,  IENEM,  HOT,  HATD,  HA3M,  IIAIS 

INTEGER  ILCND,  HCNM,  ILCNS 

INTEGER  D_IAT,D_LONG,LATS,ICNS,OBSLDC(8) 

INTEGER  EL_SIZ ,  EL_REC ,  IM_SIZ ,  IM_REC 
INTEGER  ILAD,  HAM,  HAS ,  HDD,  HCM,  HDS 
INTEGER  FLAD,  FLAM,  FLAS,FL0D,FLCM,  FLOS 
C 

REAL*4  I_VECIOR 

REAL*8  X1,Y1,Z1,X2,Y2,Z2,M(3,3) 
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REAL*8  FOCUS, XFRAME,YFRAME 
REAL*8  HEIGHT, FWR, SCALE 

set  up  graphics  environment 

create  a  color  map  with  25D  entries 

DATA  VCMJSIZE/251/ 

VOI_IIWIIS$(3*EA!IEjOOLOR_MAP  (VCM_SIZE) 

make  a  virtual  display  16.8  an  x  16.8  an 

VD_IIMJIS$CREAIE_DISFIAY{  0.0, 0.0, 512. 0,512.0, 16.8,16.8,VCM_JCD) 

fill  in  the  color  map 

DO  50  1=0,250 
I_VECT0R=I/250 . 

CALL  UIS$SET_INrENSm  (VD_ID,  I ,  I_VECIOR) 

)  CONTINUE 

set  the  background  color  and  fill  pattern 

CALL  UIS$SET_INrENSriY  (VD_ID,  0,0.0) 

CALL  UIS$SET_F0NT  (VD_ID,  1,1,'  UIS$FILL_PATTERNS ' ) 

CALL  UIS$SET_FILL_PATTERN  (VD_ID,  1,1,  PATT$C_P0REGRCUND) 

start  main  routine,  read  input 

CALL  INH7T(ELFHE,EL_SIZ,EL_KEC, 

IMFTLE,  IMJSIZ ,  IMJREC,  IIAD,  ILAM,  ILAS ,  IICD,  IICM,  ILDS , 

FLAD,  FLAM,  ETAS ,  FIDO,  FILM,  FICS ,  TTTEE) 


EL_PEC=EL_REC+1 

IM__SIZ=IM_SIZ/4 

cpen  the  input  data  files 

OPEN (UNIT=1 , FILE=ELFILE , STAIUS= ' OLD '  ,ACCESS=' DIRECT' , 
REOCSRDSIZE=EL_SIZ ,  MAXREC=EL_REC) 

OPEN  (UNIT=4 ,  FII£=IMFII£  ,  STA1US= '  OIE ',  AOCESS= » DIRECT’, 
RE00RDSIZE=IM_SIZ  ,MAXREC=IM_REC) 

CALL  READ_ELEV(ILATD, ILA1M, ILAIS,  ILCND, ILCNM,ILCNS, 
D_LAT,D_LCNG,  RELEV,  IENCN,IENEM) 

CALL  READ_IMAGE  ( IMAGE ,  IM_SIZ ,  IM_KEC) 

CALL  TER_XYZ  ( ILATD ,  ILAIM,  ILATS,  ILCND,  ILONM,  ILONS, 
DJAT,  DJJONG,  RELEV,  UNEM,  IENEN,  XYZ ,  IA,  JA, 

ILAD,  ILAM,  ILAS,  ILDD,  IICM,  IIDS, 

FLAD,  FLAM,  FLAS ,  FLDD,  FLCM,  FLOS ,  IM_SIZ ,  IM_REC) 

CALL  IM_REFAVG(IA,JA,  IMAGE,  IENEM,  IENEN,  PL_AR) 
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YIRAME=29000.0 

XERAME=35000.0 

FHE_NAME='  ' 

SAVE_FIAG=0 

P0CUS=.050 

FHR=1.0 


C 

C 

60 


70 


SCAI£>=1.2 

ITOT=XENEN*IENIM 

WD_ID<JTO$CREA!IE_WINDCW (VD_ID,  ' SYS$WORKSTATICN ' ,  TITLE, 
0.0,0.0,512.0,512.0,16.8,16.8) 


GALL  OBS_IDC  (Xl ,  Y1 ,  Z1 ,  X2 ,  Y2 ,  Z2 ,  OBSLOC) 

CALL  M_QRIEN(M,X1,Y1,Z1,X2,Y2,Z2) 

CALL  NEW_U (XI, Y1,Z1,X2,Y2, Z2,TTOT,XYZ, FOCUS, M,IMAX, 
3MAY,HIDC£N,EEPIH) 

CALL  XY2U (IA, JA, IMAX, IMAY,XERAME,  YERAME, ITOT, HIDDEN) 
CALL  HIDDEN_SURF(IA,JA,IENEM,IENEN,  DEPTH,  FL_AR,VD_ID, 
.  HIDDEN,  EEPIH_KEY,WD_ID) 


C  type  the  observer  location  and  menu 
C 

80  CALL  CLRSCRN 

WRITE (6,*)  'OBSERVER  LOCATION ' 

WRETE(6,*) 

WRTIE(6,85)  'IAT:  '  ,OBSIOC(1)  ,OBSIOC(2)  ,OBSIOC(3) 
WRTTE(6,85) 'LONG:  '  ,OBSIOC(4)  ,OBSIOC(5)  ,OBSIOC(6) 
WRITE(6,90) 'HEIGHT:  ' ,OBSIOC(7) , '  HEADING:  *,OBSIOC(8) 
WRITE(6,91)  'MAGNIFICATION:  '  ,FWR,  ’  YSCALE:  ' , SCALE 
85  FORMAT  (A8, 14, 13, 13) 

90  F0KMAT(A9,I7,A12,I3) 

91  FORMAT  (A16,F4.1,A12,F4.1) 

WRTEE(6,*) 

WRITE(6,*)  'SELECT  ONE  OF  THE  FOLLOWING' 


100 

C 


WRITE(6,*)  ' 
WRITE(6,*)  ' 
WRITE(6,*) ' 
WRITE(6,*) ' 
WRITE(6,*) ' 

READ(5,100)  ANS 
FORMAT  (Al) 


(M) agnif ication ' 

(Y) scale  factor' 
(O)bserver  Location' 
(S)ave  Image' 

(E)xit' 


IF  (ANS  .NE.  'M'  .AND.  ANS  .NE.  'm')  GO  TO  105 
WRITE (6, *)  'ENTER  MAGNIFICATION  FACTOR' 
KEAD(5,*)FWR 

102  XFRAMEXS5000.0/IWR 

YFRAME=35000 . 0/  (WR*SCALE) 

GO  TO  70 
C 

105  IF  (ANS  .NE.  'Y' .AND.  ANS  .NE.  'y')  GO  TO  110 
WRITE (6,*)  'ENTER  Y-SCALE  FACTOR' 


1 


1 


m 
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READ(5,*)  SCALE 
GO  TO  102 
C 

110  IF  (ANS  .NE.  'O'. AND.  ANS  .NE.  'o')  GO  TO  115 
GO  TO  60 
C 

115  IF  (ANS  .NE.  'S'  .AND.  ANS  .NE.  's')  GO  TO  118 

IF  (FHE_NAME  .EQ.  '  ')  THEN 
WRITE (6,*)  'ENTER  THE  NAME  OF  THE  IMAGE  FUE' 
READ (5, 116)  F3UE_NAME 
ENDIF 

SAVE_FIAG=SAVE_FIAG+1 

CALL  IMAGE_SAVE  (WD_ID ,  FHE_NAME ,  SAVE_FLAG) 

GO  TO  80 

116  FORMAT  (A13) 

118  IF  (ANS  .NE.  'E'  .AND.  ANS  .NE.  'e')  GO  TO  80 
C 

120  IF  (SAVE_FIAG  .GT.  0)  THEN 

WRITE (11)  SAVE_FIAG 
ENDIF 
CIOSE(l) 

CIOSE(4) 

CLOSE ( 10, STATOS=' SAVE') 

CLOSE ( 11 , STAIUS= ' SAVE ' ) 

CALL  UIS$DEI£TE_DISPLAY  (VD_ID) 

CALL  CLRSCRN 
END 
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******************************************************* 

SUBROUTINE  INR7T  (ELFILE ,  EL_SIZ ,  EL_REC ,  IMFIIE , 
IM_SIZ ,  IM_KEC,  HAD,  HAM,  HAS ,  HDD,  IUM,  HOS , 
FIAD,FIAM,FIAS,FLDD,FLCM,FIOS,TTTIE) 


INR7TS 


ELFILE  elevation  data  file  name 
ELJSIZ  #  column  of  ELFIIE 

EL_RE C  #  rows  of  ELFILE 

IMFIIE  image  data  file  name 
IM_SIZ  #  columns  of  IMFIIE 

IM_KEC  #  rows  of  IMFIIE 

HAD,  HAM,  HAS 

HOD,  H£M,rLOS  Northwest  comer  of  image  file 
FLAD,FIAM,FIAS 

FLOD,  FILM,  FLOS  Southeast  comer  of  image  file 
NAME  file  name  for  the  input  data  file 
TITLE  header  name  for  the  image 


OOTKJTS  =  NONE 


******************************************************** 


IMPLICIT  INTEGER  (A-Z) 

CHARACTER  EIEIIE*13,IMFIIE*13,TTn£*50,ANS*l,NAME*13 

INTEGER  EL_SIZ , EL_REC,  IM_SIZ ,  IM_REC 
INTEGER  HAD, IIAM, HAS, FLAD, FLAM, FIAS 
INTEGER  HOD,  1104,  ILOS,FLOD,FLCM,  FLOS 

start  routine  here 

CALL  CLRSCRN 

WRITE(6,*) 'Do  you  wish  to  create  an  input  data  file  (Y/N) ?  ' 
READ(5, 100)  ANS 
30  FORMAT  (Al) 

IF  (ANS  .EQ.  'Y'  .OR.  ANS  .EQ.  'y')  GO  TO  200 
IF  (ANS  .EQ.  'N'  .OR.  AN S  .EQ.  'n')  GO  TO  500 
GO  TO  5 

create  input  file 
30  CALL  CLRSCRN 

WRITE(6,*) 'Input  the  filename  for  the  input  file(*.dat):  ' 
READ(5,205)  NAME 

OPEN  (UNIT=1,NAME=NAME,STA!IUS='NEW' , 

AOCESS=' SEQUENTIAL'  ,FORM='FOFMA1TED' ) 


get  elevation  file  data 
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WRITE(6,*)  ’Enter  the  elevation  data  file  name:  ' 
READ(5,205)  ELFTLE 

WRITE(6,*) 'Elevation  file  number  of  columns:  ' 
READ(5,*)  EL_SIZ 

WRITE(6,*) 'Elevation  file  number  of  rows:  ' 

KEAD(5, *)  EL_REC 

get  image  file  data 

WRITE (6,*)  ’Enter  the  image  data  file  name:  ' 
KEAD(5,205)  IMFILE 

WRITE(6,*) 'Image  file  nurrber  of  columns:  ' 

READ(5,*)  IM_SIZ 

WRITE(6,*) 'Inage  file  number  of  rows:  ' 

READ(5, *)  IM_REC 

northwest  corner  of  image  latitude  and  longitude 

WRITE (6,*) 'Enter  the  North-west  lat/long  of  image' 
WRITE (6, *)  'Input  the  latitude  (DEC, MEN, SEC) :  ' 
READ(5,  *)  IIAD,  ILAM,  IIAS 

WRITE (6,*)  'Input  the  longitude  (DEG, MEN, SEC) :  ' 

READ  (5,  *)  ILOD,  ILCM,  ILDS 

southeast  comer  of  image  latitude  and  longitude 

WRITE (6,*)  'Enter  the  South-East  lat/long  of  image' 
WRITE (6,*) 'Input  the  latitude  (DEG, MIN, SEC) :  ' 
READ(5, *)  FTAD, FLAM, FLAS 

WRITE (6,*) 'Input  the  longitude  (DEG, MIN, SEC) :  ' 
READ(5, *)  FLOD,FLEM,FIOS 

WRITE(6, *) 'Input  a  title  for  the  image  (50  char  max. 
EEAD(5,206)  TITLE 

write  data  to  data  file 

WRITE(1,205)  ELFTLE 
WRITE (1,207)  EL_SIZ 
WRITE (1,207)  EL_REC 
WRITE (1,205)  IMFILE 
WRITE (1,207)  IM_SIZ 
WRITE (1,207)  IMJREC 
WRITE(1,209)  ILAD, ILAM, ILAS 
WRITE(1, 210)  ILOD, ILCM, ILDS 
WRITE (1,209)  FLAD, FLAM, FLAS 
WRITE (1,210)  FIDO,  FLCM,  FLOS 
WRITE(1,206)  TITLE 
C 

205  FORMAT  (A13) 

206  FORMAT  (A50) 

207  FORMAT  (13) 

2 09  FORMAT(3I3) 


nno 


210  PORTS’ (14 ,213) 

GO  TO  600 

read  an  input  file 

500  CALL  CLRSCRN 

WRITE (6, *)  'Enter  the  name  of  the  input  data  file: 
READ(5,205)  NAME 

OPEN (XJNIT=l,NAME=4lAME,SriAIUS=' OIL' , 

AOCESS=' SEQUENTIAL'  ,  FOFW=' FORMATTED' ) 
READ(1,205)  ELFILE 
READ (1,207)  EL_SIZ 
READ (1, 207)  EL_REC 
READ (1,205)  IMFILE 
READ (1,207)  IM_SIZ 
READ (1,207)  IM_REC 
READ (1,209)  HAD,  ILAM,  ILAS 
READ(1, 210)  HDD,  ILCM,  ILOS 
READ (1,209)  FIAD,FLAM,FIAS 
READ(1,210)  FIDO, FIDM, FLOS 
READ(1,206)  TITLE 
C 

600  CLOSE  (UNIT=1) 

RETURN 

END 


o  o  o 


1 


i 


******************************************************************* 


SUBROUTINE  READ_EI£V(IIAro,IIAIM>IIJ^>IIi^,IIi^,IIi»JS,  1 

D_LAT,  D_LCNG,REIEV,  IENEN,  IENEM) 

C 

C  HUS  SUBROUTINE  READS  THE  ELEVATION  FILE  AND  PLACES  THE  DATA 

C  INTO  A  REAL  50X50  ARRAY  RELEV () . 

C  INPUTS  =  NONE 

C  I 


c 

OUTPUT  =  ILATD 

southwest  latitude  of  elev  data  (degrees) 

c 

HAIM 

southwest  latitude  of  elev  data  (minutes) 

c 

HATS 

southwest  latitude  of  elev  data  (seconds) 

c 

ILOND 

southwest  longitude  of  elev  data  (degrees) 

c 

ILCNM 

southwest  longitude  of  elev  data  (minutes) 

c 

HONS 

southwest  longitude  of  elev  data  (seconds) 

c 

D  IAT 

grid  size  in  minutes  of  latitude 

c 

D  LONG 

grid  size  in  minutes  of  longitude 

c 

RELEV  () 

elev  data 

c 

IENEM 

#  rows  of  elev  data 

c 

IENEN 

#  columns  of  elev  data 

C 

c  ******************************************************************* 

IMPLICIT  INTEGER  (A-Z) 

C 

COMMON  /EIEV/ROW_EIEV/OOL_EI£V 
INTEGER  IENEN,  IENEM,  D_IAT,  D_LONG 
INTEGER  IIATD, IIA3M, ILATS,  ILOND, IIONM, IIONS 
C 

REAL  RELEV  (ROW_ELEV,aOL_EI£V) 

C 

READ(1,REC=1)  ILATD,  ILA1M,  ILATS ,  ILOND,  IICNM,  IIONS, 

D_IAT,  D_IONG,  IENEN,  IENDM 
C 

DO  10  M=l, IENEM 

READ(1,REC=M+1)  (RELEV (M,N)  ,N=1,HNEN) 

10  CONTINUE 

RETURN 
END 
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************************************************************ 

SUBROUTINE  READ_IMAGE ( IMAGE , IM_SIZ , IM_REC) 

HUS  SUBROUTINE  READS  THE  IMAGE  DATA  INTO  AN  ARRAY. 
INPUTS  =  IM_SIZ  number  of  columns  of  image 
IM_REC  number  of  rows  of  image 

OUTIUr  =  IMAGE  ()  image  gray  level  file 

C************************************************************ 

IMPLICIT  INTEGER  (A-Z) 

CCflMDN  /IMAGE/ROW_IMAGE,OOL_IMAGE 
C 

BYTE  IMAGE (ROW_IMAGE,OOL_IMAGE) 

C 

INTEGER  IM_SIZ , IM_REC 
C 

DO  10  IR=1,IM_REC 

READ(4,REC=IR)  (3MAGE(IR,IC)  ,IC=1,IM_SIZ*4) 

10  CONTINUE 
REIURN 
END 
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*************************************************************** 

SUBROUTINE  TER_XYZ(IAID,LAIM,  LATS,  IIOND,  LONM,  LONS, 

.  D_LAT,D_LCNG,REIL£V,IENDM,IENCN,XYZ,IA,JA, 

.  ILAD,IIAM,IIAS,ILOD,ILCM,IIOS, 

.  FLAD,  FLAM,  FLAS ,  FIDO,  FLCM,  FIDS ,  IM_SIZ ,  IM_REC) 

THIS  SUBROUTINE  CONVERTS  EACH  ELEVATION  POINT  TO  ITS  DMS 
EQUIVALENT.  IT  USES  THE  FACT  THAT  EACH  POINT  REFRESEMS 
AN  EQUAL  CHANGE  FRCM  THE  LAST.  THE  REFERENCE  IAT/IONG 
AND  THE  DELIA  IAT/IONG  ARE  PASSED  IN  THE  SUBROUTINE  CALL. 

. .  .ASSUMES  COLOR  COMPLETELY  OVERLAPS  THE  ELEV  DATA 

. .  .NEED  TO  KNOW  THE  IAT/LONG  OF  THE  COLOR  IMAGE 

. .  .ALSO  CALCS  THE  IA() ,  JAQ  COORDINATES  FOR  EACH  ELEV  FT 

ELEVATION  DATA 

INPUTS  LAID  -  REFERENCE  LATITUDE  (DEGREES) 

LA3M  -  REFERENCE  LATITUDE  (MINUTES) 

LATS  -  REFERENCE  LATITUDE  (SECONDS) 

IIOND  -  REFERENCE  LONGITUDE  (DEGREES) 

IONM  -  REFERENCE  LONGITUDE  (MINUTES) 

IONS  -  REFERENCE  LONGITUDE  (SECONDS) 

D_LAT  -  DISTANCE  BETWEEN  ROWS 
D_IONG  -  DISTANCE  BETWEEN  COLUMNS 
RELEV  ()  -  ELEVATION  DATA  FILE 
IENCM  -  NUMBER  OF  ROWS 
UNDN  -  NUMBER  OF  COLUMNS 

OUTPUTS  XYZ  -  OUTFUT  DATA  FILE 

*  *  *  *  *  *  *  *  *  *  ★  *  *  it  ★  ★  it  ★  *  *  *  *  *  ★  ★  *  ★  *  ★  *  ★  *  *  ★  *  *  **  irk  ★★  irk  *  ★★  *  *  *  **  icirk  * 

IMPLICIT  INTEGER  (A-Z) 

COMMON  /ELEV/  ROW_ELEV ,  OOL_ELEV ,  TOTJELEV 

INTEGER  IA  (TOT_ELEV) ,  JA  (TOT_ELEV) 

INTEGER  LATD,  IAIM,  LAIS ,  IIOND,  IIONM , D_LAT ,  D_IONG 
INTEGER  LCND,  LONM,  LONS 
INTEGER  HAD,HAM,ILAS,ILOD,IIOM,IIOS 
INTEGER  FTAD,  FLAM,  FLAS ,  FIOD,  FLCM,  FLOS 
INTEGER  3M_SIZ,IM_REC 

PEAL  RELEV  (RCW_ELEV,  OOL_ELEV) 

REAL*3  HEIGHT,  X,  Y,Z,XYZ  (1UT_ELEV,  3) 

convert  to  seconds 

CALL  CDNV2SEC(ILAD,  ILAM,  ILAS) 

CALL  OONV2SEC(ILOD,  IICM,  ILOS) 

CALL  OONV2SEC(FLAD, FLAM, FLAS) 

CALL  OONV2SEC(FLOD, FIOM,  FLOS) 
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c 


CELLAT=HAS-FIAS 

nKIJCN=FD36-IIClS 


set  flag  for  proper  hemisphere 

IF  (HAD  .GE.  0)  MEN 
IAT_FIAG=1 

VTKV. 

IAT_FIAG=-1 

ENDIF 

IF  (HDD  .GE.  0)  THEN 
IDNJF1AG=1 
FTCT! 

L0N_FIAG=-1 

ENDIF 

find  initial  seconds  of  elev.  data 

CALL  CCNV2SEC(IATO,IA!IM,  IATS) 
CALL  G3tft/2SEC(ILCND,IDNM,LONS) 

set  up  for  1st  increment 


LATS=IATS-D_IAT 
II£NS=LDMS-D  LONG 


start  routine  here 
IR=0 

DO  10  M=i,imn 
IAIS=1ATS+D_IAT 
IONS=IIONS 
DO  20  N=1,IENEN 
IR=IRfl 

ION^IONS+D_IDNG 
HEHGHT=RELEV  (M,N) 

BCM=mr(  (3M_REC-1)  *  (IIAS-LAIS) /DELLAT)  +1 
001>INr(  ( (IM_SIZ*4)  -1)  *  (LONS-ILQS) /DETION)  +1 

the  rcw  and  column  coordinates  are  1-512 

IA(IR)=OOL 
JA(IR)=R0W 

CALL  EMS2XYZ  (0,0, IAIS,0,0,IDNS,HEiair,X,Y,Z) 
XYZ(IR,1)=X 
XYZ(IR,2)=Y 
XYZ(IR,3)=Z 
20  CONTINUE 

10  CONTINUE 
REIURN 
END 
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******************************************************************** 
SUBROUTINE  CM52XYZ  (IATO,IA3M,IATS,  K*TO,LCNM,  LONS,  HEIGHT,  X,Y,  Z) 


THIS  SUBROUTINE  CONVERTS  CMS  DATA  TO  X,Y,Z. 


INRJTS  =  LAID 

latitude  in  degrees 

IA3M 

latitude  in  minutes 

LATS 

latitude  in  seconds 

KM) 

longitude  in  degrees 

KNM 

longitude  in  minutes 

KNS 

longitude,  in  seconds 

HEIGHT 

height  above  sea  level  (meters) 

OUTH7T  =  X,  Y,  Z 

XYZ  coordinates  of  CMS  point 

IMPLICIT  INTEGER  (A-Z) 

C 

INTEGER  IATO,LAlM,LATO,KND,LONM,LONS 
C 

REAL*S  PHI,  IAMDA,N,X,Y,Z,  HEIGHT 
REAL*8  PI ,  Cl ,  C2 ,  C3 ,  E_SQUARE ,  A ,  RADIAN 
C 

PARAMETER ( PI=3 . 14159265358793238) 

PARAMETER  (Cl=180 .  ,C2=60.  ,C3=3600. ) 
PARAMETER  (E_SQUAHE=0 . 006768658 ,  A=6378206 . 4 ) 
C 

RADIAN=PI/C1 

C 

IF  (LAID  .GE.  0)  THEN 
PHI=  (IAIDf  IAUVC2+IATS/C3 )  *RADIAN 
ELSE 

PHI=  (IATD-LATM/C2-LAIS/C3 )  *RADIAN 
END  IF 
C 

IF  (KM)  .GE.  0)  THEN 
LAMDAf  (KM>KMVC2+ICNS/C3 )  *RADIAN 
FT.9F. 

LAMDA=(KM>K3N^C2-KXIS/C3)  *RADIAN 
END  IF 
C 

N=A/SQFT  (1-E_SQUARE*SIN  (PHI)  *SIN  (HU) ) 
X=(N4-HEK2fT)  *C0S  (FKE)  *00S  (LAMDA) 

Y=  (N+HEIGHT)  *CCS  (FKC)  *SIN  (IAMDA) 

Z=  (N* ( 1-E_SQUARE)  +HEIGHT)  *SIN(FHI) 

RETURN 

END 
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c 

0************************************************************* 

c 

SUBROUTINE  IM_KEFAVG (IA,  JA,  IMAGE,  IENEM,  IENEN, PL_AR) 

HUS  SUBROUTINE  OCNSTKJCTS  THE  GECMEIKY  FH£ 

THAT  CONTAINS  THE  THREE  NODAL  POINTS  THAT  MAKE  UP 
AN  IMAGE  PLANE  AND  THE  GREY  SCALE  VAHJE  ASSIGNED  TO 
THAT  PLANE. 

INPUTS  =  IMAGE  ()  image  gray  level  file 

IENEM  #  rows  of  elev  data 

IENEN  #  columns  of  elev  data 

IA()  image  column  index 

JA()  image  row  index 

OUnvr=  PL_AR()  array  holding  nodes  and 

gray  value  for  each  panel 

*  *  -k  irk  ★  *  *  *  ★  *  ★  it  *  it  ★  *  *  ★  *  *  *  *  *  ★  *  ★  *  *  *  it  -k  ★  *  *  *  *  *  *  *  *  *  *  *  ★  ★  it  *  *  ★  ★  ★  ★  it  it  *  *  ★  * 

IMPLICIT  INTEGER  (A-Z) 

CCMM3N  /EI£V/RCW_EIEV,  00L_ELEV,1UT_ELEV,NPIANES 
OOMCN  /IMAGEyROW_IMAGE ,  OOL_IMAGE 

BYTE  IMAGE  (FCW_IMAGE ,  OOL_IMAGE) 

INTE3GER  IA  (TUTJELEV)  ,  JA  (TOT_ELEV)  ,  PLAR  (NFLANES  r-  5) 

INTEGER  IY,M,N,  IGREY1 ,  IGREY2 ,  L,  LI 
INTEGER  IENEM,  IENEN, NODE_A,NODE_B,NODE_C,NODEJD,  IR 
INTEGER  IOCUNri,IOOUNI2,irOTl,ITOr2 

REAL*8  SLOPE, YINT 

DO  90  IR=1 ,  HNEM-1 
ID=(IR-1)  *IENEN 
DO  80  N=1+IL,  IENEN-l+IL 
NODE_A=N 
NODE_B=N+l 
NOEE_C=N+IENEN 
N0DE_D=N0DE_C+1 

SIOFE=  ( JA (N0DE_C)  -JA(N0DE_B) )  *1.0/  (IA(NODE_C)  -IA(NODE_B) )  *1. 0 
YINr=l .  0*JA  (N0DE_B)  -SLOPE*  IA  (NOCE_B) 

new  knew  the  slope  and  y-intecept  of  the  line  forming  the  triangle 

nun=o 
nOT2=0 

I00UNTl=O 
I0dUNT2=0 

DO  70  M=IA(NODE_B)  ,IA(NODE_A)  ,-l 
IY=SLOPE*MfYINr 
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DO  50  L=JA(NODE_B)  ,IY,-1 
IF  (IMAGE  (L,M)  .Iff.  0)  THEN 
nOTl=rron+IMAGE  (L,M)  +256 

iron=iTan+iMAGE  (l,m) 

END  IF 

I(XUNTl=ICOUNn+l 

CONTINUE 

3F(M.Iff.IA(N0QE_B)  )THEN 
DO  60  I?=IY ,  JA  (NOQ£_C )  ,-l 
IF  (IMAGE(L,M)  .Iff.  0)  THEN 
ITOT2=ITOT2+IMAGE(L/M)+256 
FT  SR 

rror2=iTor2+iMAGE  (l,m) 

END  IF 

I<XtJNT2=IC0UNT2+l 
CONTINUE 
END  IF 
CONTINUE 
LI= (N- ( IR-1 ) ) *2 
I<3REXl=ITOri/ICaUNTl 
IGEEY2=ITOT2  /  I00UNT2 
PL_AR(L1-1 , 1)  =NOCE_A 
PL_AR(L1-1 , 2)=N0DE_B 
PLAR  (Ll-1,3)  =NOCE_C 
FL_AR(L1-1, 4) =IGREY1 
PL  AR(U,1)=N0QE_B 
PL_AR(L1,2)=N0DE_D 
PL_AR(L1,3)==N0CE_C 
PLAR  (LI,  4)  =IGREY2 
OCNTINUE 
OCNTINUE 
RETURN 


vo  oo  (-><■>  oooo  ooooooooooooooo 


************************************************************ 

SUBROUTINE  0BS_L0C(X1,Y1,Z1,X2, Y2,Z2,OBSIDC) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  OBSERVER  XYZ 
LOCATION  FROM  DESIRED  LAT.  AND  LONG.  INPUTS. 

INPUTS  =  NONE 

OUTFUIS  =  XI, Yl, Z1  observer  location 

X2,Y2,Z2  observer  location  direction 

OBSIOCO  observer  location  array 

************************************************************ 
IMPLICIT  INTEGER  (A-Z) 

CHARACTER  ANS*1 

INTEGER  IATD,  IAIM,  LAIS ,  LOND,  LCNM,  IONS ,  CTS 
INTEGER  OBSLOC(8) 

REAL*8  XI, Y1,Z1,X2,Y2,Z2, HEIGHT 

read  in  last  values  of  observer  location 

IATD=OBSIOC  (1) 

IATM=OBSIOC(2) 

LATS=OBSIOC  ( 3 ) 

LOND=OBSIOC(4) 

ICNM=OBSIOC  (5) 

LONS=OBSLOC(6) 

HEIGHT=OBSIOC  (  7 ) 

CTS=OBSIOC(8) 

FIAG=0 

CALL  CLRSCRN 

WRITE (6,*)  'OBSERVER  LOCATION' 

WRITE(6,*) 

WRTTE(6,8) 'LAT:  '  ,OBSIOC(1)  ,OBSIOC(2)  ,OBSLOC(3) 
WRITE(6,8) 'IONG:  ' ,OBSIOC(4) ,OBSIOC(5) ,OBSLOC(6) 

WRITE (6, 9) 'HEIGHT:  ' ,OBSIOC(7) , '  HEADING:  \OBSIOC(8) 
FORMAT  (A8, 14, 13, 13) 

FORMAT  (A9, 17,  A12, 13) 

WRITE(6,*) 

IF  (FLAG  .EQ.  1)  THEN 
WRITE(6,*)  ARE  THERE  ANY  CORRECTIONS' 

WRTTE(6,*) 

ENDIF 

WRITE (6,*)  'SELECT  ONE  OF  THE  FOLLOWING' 

WRITE{6,*)'  (1)  latitude' 

WRITE(6,*)'  (2)  Longitude' 

WRITE (6,*)'  (3)  Height ' 


WRITE(6,*)'  (4)  Heading* 

WRITE(6,*)'  (5)  All* 

n 

WRITER,*)*  (6)  None* 

READ(5,100)  ANS 

100 

PCfFMftT(Al) 

n 

IF  (ANS  .BQ.  *6')  GO  TO  160 

IF  (ANS  .NE.  *1*  .AND.  ANS  .NE.  '5') 

GO  TO  60 

WRITE  (6,*)  'INPUT  OBSERVER  IATI'iUDE 

-DEGREES: 

READ(5,*)IAID 

WRITE(6,*) ' 

-MINUTES: 

READ(5,*)IA1M 

WRITE(6,*)  * 

-SECONDS: 

c 

READ(5,*)IATS 

60 

IF  (ANS  .NE.  '2'. AND.  ANS  .NE.  *5') 

GO 

TO  70 

WRITE(6, *) 'INPUT  OBSERVER  LONGITUDE 

-DEGREES: 

READ(5,*)IDND 

WRITE(6, *) ' 

-MINUTES: 

READ(5,*)IDNM 

WRITE(6,*) ' 

-SECONDS: 

c 

READ(5,*)IDNS 

70 

IF  (ANS  .NE.  *3 ' .AND.  ANS  .NE.  *5') 

GO 

TO  80 

WRITE (6,*)  'INPUT  OBSERVER  HEIGHT 

-METERS:  ' 

r 

READ  (5,*)  HEIGHT 

w 

80 

IF  (ANS  .NE.  *4'  .AND.  ANS  .NE.  *5') 

GO  TO  105 

WRITE(6, *)  'INPUT  THE  COURSE 
READ(5,*)CTS 

DEGREES: 

105 

IF  (ANS  .IT.  .OR.  ANS  ,GT.  '6') 

GO  TO  1 

C  load  the  values  into  the  obs_loc  array 
C 

150  OBSLDC  ( 1 )  =IATD 

OBSLDC ( 2 ) =IATM 
GBSLOC ( 3 ) =LATS 
OBSLDC (4 )=LOND 
OBSIDC  (5)  =IDNM 
OBSIDC  ( 6)  =LONS 
OBSIDC  (7)  =HEIGHr 
OBSLDC  (8)  =CIS 
FLAG=1 
GO  TO  1 
C 

C  convert  to  XYZ  coordinates 
C 

160  CALL  CMS2XYZ  (LAID,  LAlM,IA!rS,IOND,LONM,  IONS,  HEIGHT,  XI,  Y1,Z1) 
C 

C  calc  new  position  based  on  course 
C 
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IF  (LAID  .GT.  0)  THEN 

IATS=5*OOS  (2*3 . 14159*CIS/360)  +IA3S 
FT  .OF. 

IATS=IATS-5*SIN  (2*3 . 14159*CIS/360) 

END  IF 

IF  (LOND  .GT.  0)  THEN 

LCNS=5*SIN(2*3 . 14159*CIS/360)  +LCNS 
FTfiF 

I£NS=L0NS-5*SIN  (2*3 . 14159*CTS/360) 

END  IF 

convert  new  position  to  XYZ  coordinates 

CALL  EMS2XYZ  (IATD,  IA3M,  LAIS ,  I£ND,  LCNM,  ILNS ,  HEZQfT,  X2 ,  Y2 ,  Z2 ) 

RETURN 

END 
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c  ************************************************************ 
SUBROUTINE  M_CRIEN(MfXl,Yl,Zl,X2,Y2,Z2) 

THIS  ROUTINE  EETEFMINES  TOE  ANGLES  OF  ROTATION 
OF  THE  IMAGE  PLANE  WITH  RESPECT  TO  THE  WORLD  COORDINATE  AXEES. 
AND  CALCULATES  THE  NEW  M  MATRIX  FOR  THE  NEW  VIEWER  LOCATION 

INFUTS  =  XL,Y1,21  observer  location 

X2,Y2,Z2  observer  direction  point 

OUTFUIS  =  M()  3D  -  2D  xform  matrix 

*************************************************************** 
IMPLICIT  INTEGER  (A-Z) 

KEAL*8  MAGN_X,MAGN_Y,MAGN_Z 

REAL*8  X_VECX ,  X_VECY , X_VECZ , Y_VECX , Y_VECY ,  Y_VECZ 
REAL*8  Z_VECX,Z_VECY,Z_VECZ,X1,Y1,Z1,X2,Y2,Z2 
REAL*8  M(3,3) 

get  the  two  OBS  IOC  points 
and  set  up  vectors 

Y_VECX=X1 
Y_VECY=Y1 
Y_VECZ=Z1 

select  another  point  along  the  track 

Z_VECX=X1-X2 
Z_VECY=Y1-Y2 
Z_VECZ=Z1-Z2 

use  the  cross  product  of  Y  CROSS  Z  to  obtain  the  X  vector 

X_VECX=(  ( Y_VECY*Z_VECZ )  -  (Y_VECZ*Z_VECY) ) 

X_VECY=(  ( Y_VECZ*Z_VECX)  -  (Y_VECX*Z_VECZ ) ) 

X_VECZ=(  (Y_VECX*Z_VECY)  -(Y_VECY*Z_VECX) ) 

MAGN_Z=SQFT(  (ZJTXX**2)  +  (ZJ®CY**2)  +  (Z_VECZ**2) ) 

MAGN_X=SQRT ( (X_VECX**2 )  +  (X_VECY**2 )  +  (X_VECZ**2 ) ) 

MAGN_Y=SQRT(  (YJ/ECX**2)  +  (Y_VECY**2)  +  (Y_VECZ**2) ) 

M(l, 1) =XJ/ECX/MAGN_X 
M  ( 1 , 2 )  =X_VECY/MAQ^_X 
M(l, 3) =X_VECZ/MAGN_X 
M(2 , 1)  =Y__VECX/MAGN_Y 
M(2 , 2)  =Y_VECY/MACN_Y 
M(2, 3)=Y_VECZ/MAGN_Y 
M(3, 1)  ^_VECX/MAGN_Z 
M(3 , 2)=Z_VECY/MAQJ_Z 
M  ( 3 , 3 )  =Z_VECZ/MAGN_Z 
PERM-' 
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**************************************************************** 

SUBROUTINE  NEW_U(X1,Y1,Z1,X2, Y2,Z2(ITOT, 

XYZ ,  FOCUS  ,M,  IMAX,  IMAY ,  HIDDEN,  DEPIH) 

CHS  SUBROUTINE  CCMFUIES  THE  NEW  IA  AND  JA  SCREEN 
COORDINATES  FROM  THE  GP/EN  OBSERVER  IOCATICN. 

INPUTS  =  X1,Y1,Z1  observer  location 

X2,Y2,Z2  observer  direction  point 

ITOT  #  elev  points 

M()  3d  -  2d  xform  matrix 

FOCUS  focal  length 

XYZ()  dms  coordinates  of  elev  data 

CUTFUTS  =  IMAX()  x  image  coordinate 

IMAYQ  y  image  coordinate 

HIDDEN  ()  flag  for  points  hidden  behind  FOV 
DEPIHQ  depth  of  each  point 

*************************************************************** 
IMPLICIT  INTEGER  (A-Z) 

CCMM3N  /ELEV/  ROWJSEEV,  OOL_ELEV,  TUTJEXEV 

INTEGER  ITUT,  LR, HIDDEN  (TOT_EI£V) 

REAL*8  X1,Y1,Z1,FOOUS,M(3,3) 

REAL*8  X,Y,Z,XIMA,Y3MA,DENCM,DEPm(TOT_ELEV) 

REAL*  8  XYZ  (TOT_ELEV,  3 )  ,  IMAX  ('IUT_EIEV) ,  IMAY  (TOT_ELEV) 

REAL*8  X2,Y2,Z2, OBSJVECX , OBS_VECY , OBS_VECZ 
REAL*8  OBJ_VECX ,  OBJVECY ,  OBJ_VECZ ,  MAG_OBS ,  MAG_OBJ 
REAL* 8  OC6JIHETA, VALUE 

VAIDE=90.0 

VAIUE=COS  (VALUE/ 57 . 29577951) 

TO  20  IR=l,nOT 
X=XYZ(LR,1) 

Y=XYZ(IR,2) 

Z=XYZ (IR, 3) 

calc  the  cbs_vec 

OBS_VECX=X2-Xl 
OBS_VECY=Y2-Yl 
OBS_VECZ=Z  2 -Z 1 

MAG_OBS=SQRT  (OBS_VECX**2+OBS_VECY**2+OBS_VBCZ**2) 
calc  the  cbj_vec 
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0BJ_VECX=X-X2 
0BJ_VECY=Y-Y2 
0BJ_VECZ=Z-Z2 

MAG_OBJ=SQRT  (OBJJ/ECX**2^BJ_VECY*^2+OBJ_VECZ**2 ) 

calc  cos  theta,  where  theta  is  angle  between  line  of  sight 
and  object  point. 

OX>_THETAOBS_VECX*OBJ_VECX+OBS_VECY*OBJ_VECY+OBS_VECZ*OBJ_VECZ 
OOSJIHEIA<X3S_lHEnV(^_OBS*MftG_OBJ) 

now  check  if  CDS  (theta)  >  0  =>  angle  <  90  degrees 

IF  (CDSJIHEIA  .GT.  VAKJE)  THEN 

DENCH=M(3 , 1)  *  (X-X1)4M(3 , 2)  *  (Y-Yl)  4M(3 , 3)  *  (Z-Zl) 
X]MAfHFOCUS*(M(1,  1) * (X-Xl) 4M(1, 2) * (Y-Y1)+M(1, 3)  *(Z-Z1)  )/EENCM 
YIMft=-RXUS*  (M(2 , 1)  *  (X-Xl)  -HM(2 ,  2)  *  (Y-Yl)  -H«(2 ,  3)  *  (Z-Zl)  )/DEN0M 
X3MA=XIMA* 1000000 . 0 
YIMA=YTMA*1000000 . 0 

save  xima,  yima  in  IMAX() ,  IMAY()  array  temporarily 

IMAX  (IR)  =XIMA 
IMAY ( IR) =YIMA 

EEPIH(IR)=SQRT(((X1-X)  **2)+( (Yl-Y) **2)+( (Zl-Z) **2) ) 

HIDDEN (IR)=0 
EISE 

HIDDEN  (IR)=1 
END  IF 
20  CONTINUE 
RETURN 
END 
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***************************************************************** 


SUBROUTINE  XY2U(IA,JA, 

IMAX ,  IMAY ,  XFRAME ,  YFRAME ,  ITOT ,  HIDDEN) 

THIS  SUBROUTINE  TAKES  THE  IMAGE  POINTS  XIMA,  YIMA  AND 
OCNVERIS  THEM  TO  I,J  SCREEN  POINTS 
THIS  IS  THE  AFFIN  XPORM  —  SEE  REF  3 


INPUTS  =  IMAX() 

x  image  coordinate 

3MAY() 

y  image  coordinate 

XFRAME 

image  coordinate  x  axis  size 

YFRAME 

image  coordinate  y  axis  size 

rror 

number  of  elev  points 

HIDDEN  () 

hidden  point  flag 

OUTPUTS  =  IA() 

screen  x  coordinate 

JA() 

screen  y  coordinate 

****************************************************************** 
IMPLICIT  INTEGER  (A-Z) 

OCTM3N  /ELEV/  RCW_ELEV ,  CDL_EI£V ,  TOT_ELEV 

INTEGER  IA(TUT_ELEV) ,  JA(TOTJSXEV) ,  HIDDEN  (TOT_ELEV) ,  TIUT 

KEAL*8  XIMA,  YIMA,A1,  A2  ,B1,B2,C1,C2,D£NCM 
REAL*8  IMAX  (TOT_ELEV) ,  IMAY  (TUT_EI£V)  ,  XFRAME ,  YFRAME 

DATA  I  MAX, J_MAV512 .0,512.0/ 

calc  the  affin  parameters 

A1=XFRAME/ ( I_MAX*1 . 0) 

A2=0.0 
B1=0.0 

B2=YFRAMEy  ( J_MAX*1 . 0) 

Cl=-XFRAME/2 . 0 
C2=-YFRAME/2.0 

DO  25  IR=1,ITOT 

IF  (HIDDEN (IR)  .EQ.  1)  GO  TO  25 
XIMAp=IMAX(IR) 

YIMA=IMAY ( IR) 

IA(IR)  =  (XIMA-C1)/A1 
JA(IR)  =  (YIMA-C2)/B2 
25  OCNTINUE 

REIURN 
END 
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SUBROUTINE  HIDCEN_SURF(IA,JA,IENEM/IENEN,EEPIH, 

PL_AR,  VD_ID ,  HIDCEN ,  EEPIH_KEY,WD_ID) 

CHS  ROUTINE  WILL  CHECK  EACH  PANEL  AGAINST  AIL 
OTHERS  FOR  OVERLAP. 

INPUTS  «=>  IA(  ) 

JA(  ) 

IENEN 

IHCH 

HIDCEN  ()  flag  for  points  hidden  behind  POV 


OUTPUTS  ■  NONE 

a*********************************************,**************** 

IMPLICIT  INTEGER  (A-Z) 

OCMCN  /ELEV/  ROW_ELEV ,  CDL_EL£V ,  TUT_ELEV ,  NPLANES 

INTEGER  K3?EY,NOCE_A,NM3E_B<NCCE_C 
INTEGER  IENEN, HNEMrIR,I,J 

INTEGER  IPLANES , HIDDEN (TDTJELEV)  , FL_AR(  NPLANES, 5) 

INIEGER  IA(TOTJELEV) ,  JA(TOT_ELEV)  ,EEPIH_KEY  (NPIANES) 

INTEGER  X1,Y1,X2, Y2,X3,Y3 
KEAL*8  DEPIH  (TOT_ELEV) 

LOGICAL  SORTED, PLUS, MINUS 

calculate  the  number  of  planes 

IPLANES=( (IENEM-1) *2) * (IENDN-1) 

C0UNP=0 

read  in  the  data  for  each  plane 

DO  100  IR=1,  IPLANES 
PTA=PL_AR(IR, 1) 

PTB=PL_AR(IR,2) 

PK>=PL_AR  ( IR,  3 ) 

IF  (HIDCeN(PEA)  .EQ.  1  .CR.  HIDDEN  (PTB)  .EQ.  1  .OR.  HIDDEN (PTC) 
.EQ.  1)  THEN 
GO  TO  100 
FTfiF 

OOUNI^KDOUNIM-l 
EEPm_KEY  (COUNT)  =IR 
ENDIF 
DO  CONTINUE 

sort  depth  values 


c 

IPtANESKXJUNT 
DO  105  LfI, IPLANES 
IR=CEPIH_KEY  (L) 

PL_AR (IR,  5)  =MAX(EEFIH (PL_AR ( IR,  1) ) , 
EEPIH(FL_AR(IR,  2) )  ,  DEPTH (PL_AR(IR,  3) ) ) 

105  CONTINUE 
C 

CALL  QUICKSCRT  (PL_AR,  IPLANES ,  EEPIH_KEY) 

C 

CALL  UISDC$ERASF  (WD_ID) 

C 

C  new  begin  selection 
C 

109  DO  110  K=IFLANES ,1,-1 

LR=CEPm_KEY(K) 

C 

C  draw  to  virtual  display 
C 

NOEEA^PLAR  ( IR ,  1) 

NOTE_B=PL_AR(IR>  2) 

NCCE_C=PL_AR(IR,  3) 

IGKEY=PL_AR  ( IR ,  4 ) 

C 

C  limit  the  gray  scale  to  250  shades  of  gray 
C  vice  256  shades 
C 

IF  (IGREY  .EQ.  0)  IGREY=1 
IF  (IGREY  .GT.  250)  IGREY=250 
X1=IA  (N0DEA) 

Y1=JA(N0DE_A) 

X2=IA  (N0DE_B) 

Y2=JA  (NODEB) 

X3=IA(N0DE_C) 

Y3=JA(NOCE_C) 

CALL  UIS$SET_WRITING_INEEX  (VD_ID,  1,1, 1®EY) 
CAUL  UISDC$PI0T(WD_ID,  1,X1, Y1,X3, Y3,X2,Y2) 
C 

110  CONTINUE 
REIURN 
END 
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********************************************** 
SUBROUTINE  CLRSCRN 
Dlls  routine  clears  the  screen 


Inputs  =  NONE 


Outputs  »  NONE 


********************************************** 
IMPLICIT  INTEGER  (A-Z) 

CALL  SMS$CKEATE_PASTEBOARD(PB_ID) 

CALL  SM3$ERASE_PASTEBQARD(PB_3D) 

RETURN 

END 


***************************************************** 

SUBROUTINE  OONV2SEC(DBG,MIN,SEC) 


INPUTS  =  DBG 
MIN 
SEC 


OUT FUTS=  SEC 

***************************************************** 
IMPLICIT  INTEGER  (A-Z) 

IF  (DBG  .GE.  0)  THEN 
SEC=SEC4MIN*60+DBG*3600 

FTfiK 

SEO=QBG* 3 600-KIN*  60-SEC 
ENDIF 
RETURN 
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Q  ************************************************** 

SUBROUTINE  QUICRSCRr (ARRAY, COUNT, KEY) 

INR7T  =  ARRAY  ()  array  to  sort,  sort  on  5th  field 
KEYQ  key  index  to  main  array 
OOUNT  number  of  elements  to  sort 

OUTEUr=  KEY()  new  key  index  to  main  array 

**************************************************** 

IMPLICIT  INTEGER  (A-Z) 

CX>M3N/ELEV/NN,I^l,'IUr,NPLANES 
INTEGER  ST (100, 2) 

INTEGER  COUNT, KEY  (NPLANES) 

INTEGER  ARRAY(NPLANES,5) 

C 

SP=1 

ST(1, 1)=1 

ST(1,2)=OOUNT 

DO  WHILE  (SP  .NE.  0) 

D=ST(SP,1) 

R=ST(SP,2) 

SP=SP-1 

DO  WHILE  (R  .CT.  L) 

LI=L 

RI=R 

3NDEX=INr(  (L*R)/2) 

SA=ARRAY (KEY  (INDEX) ,  5) 

DO  WHILE  (LI  .LE.  RI) 

DO  WHILE  (ARRAY (KEY (LI), 5)  .IT.  SA) 

LI=LI+1 
END  DO 

DO  WHILE  (ARRAY (KEY (RI), 5)  .GT.  SA) 

RI=RI-1 
END  DO 

IF  (LI  .LE.  RI)  THEN 
TEMP=KEY(LI) 

KEY (LI)=KEY (RI) 

KEY(RI)=TTMP 
U=LI+1 
RI=RI-1 
ENDIF 
END  DO 

IF  ( (R-LI)  .GT.  (RI-L) )  THEN 
IF  (L  .IT.  RI)  THEN 
SP=SF4-1 
ST(SP,  1)  =L 
ST(SP,2)=RI 
ENDIF 
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C*****w******************************************* 

c 

SUBROUTINE  IMAGE_SAVE  (WD_ID,  FILE_NAME ,  FLAG) 
C 

0************************************************* 

c 

IMPLICIT  INTEGER  (A-Z) 

C 

CHARACTER  FILE_NAME*(*)  ,NAMEl*23fNAME2*23 
C 

INTEGER  WD_ID,FIAG 
C 

DATA  RETUN/O/ 

DATA  pwronyo/ 

DATA  RHEIGHT/O/ 

DATA  BPPIXEI/O/ 


determine  the  size  of  the  inage  buffer 

CALL  UISDC$READ_1MAGE (WD_ID,  0,0, 512 , 512  ,RWIIHH,KHEIGHr,  BPPIXEL, 
ENOODED1,0) 

Rm£N=iwinm*RHEi(3fr 

open  the  external  files 

IF  (FLAG  .EQ.  1)  THEN 
C1=INEEX  ( FILE_NAME , • . ' ) 

IF  (Cl  .EQ.  0)  THEN 

C2=INDEX (FILE_NAME , '  ') 

NAME1=FILE_NAME ( 1 : C2-1) // ' .SIZ' 

NAME2=FIIE_NAME(l:C2-l)//,  .IMG' 

ELSE 

NAME1=FHE_NAME ( 1 : Cl-1 )  //  '.SIZ' 

NAME2=FH£_NAME 

ENDIF 

WRITE  ( 6 ,  * )  NAME1 ,  NAME2 

OFEN  (UNIT=11 ,  FHE=NAME1 ,  STAIUS= '  NEW ' , 

AOCESS= '  SEQUENTIAL1 F3BM= '  UNFORMATTED') 


OPEN  (UNIT^10,FILE=4®ME2,STAIUS=' NEW', 

P0RM=' UNFORMATTED' ) 

WRITE  ( 11)  FWIDIH ,  FREIGHT,  BPPIXEL,  REITEN 

allocate  virtual  memory  for  the  buffer 

STAIUS=LIB$GEr_VM  (RETLEN ,  ENOOCED) 

IF  (.NOT.  STATUS)  CALL  LrB$STOP(%VAL(STAIUS) ) 
ENDIF 

extract  and  store  private  data  in  buffer 
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C 

CALL  UISDC$READ_IMAGE  (WD_ID,  0,0,512, 512 ,  FWIDIH,  RHEIGHT ,  BPPIXEL, 
%VAL(EN00CED)  ,RETL£N) 

call  subroutine  to  write  the  contents  of  the  buffer 

CALL  BUFFERWRITE  ( %VAL (ENCODED)  fRETI£N,10) 

RETURN 
END 


SUBROUTINE  BUFFERWRITE  (BUFFER,  LENGIH,  IUN) 
********************************************* 

IMPLICIT  INTEGER  (A-Z) 

C 

BYTE  BUFFER  (LENGIH) 

C 

WRITE  (IUN)  BUFFER 
C 

REIURN 

END 
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